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ABSTRACT 

We consider (/-state Potts models coupled by their energy operators. Restricting 
our study to self-dual couplings, numerical simulations demonstrate the existence 
of non-trivial fixed points for 2 < q < 4. These fixed points were first predicted 
by perturbative renormalisation group calculations. Accurate values for the cen- 
tral charge and the multiscaling exponents of the spin and energy operators are 
calculated using a series of novel transfer matrix algorithms employing clusters 
and loops. These results compare well with those of the perturbative expansion, 
in the range of parameter values where the latter is valid. The criticality of the 
fixed-point models is independently verified by examining higher eigenvalues in the 
even sector, and by demonstrating the existence of scaling laws from Monte Carlo 
simulations. This might be a first step towards the identification of the conformal 
field theories describing the critical behaviour of this class of models. 

PACS numbers: 05.50. +q,64.60.Fr,75.10.Hk,75.40.Mg 



1 Unite Mixte de Recherche CNRS UMR 7589. 

2 Laboratoire associe aux universites Paris 6, Paris 7 et au CNRS. 



1 Introduction 



It is well-established that many spin models of statistical physics with random-valued 
nearest-neighbour couplings possess distinct critical points, their critical properties being 
different from those of the corresponding fixed-coupling models. To ensure the existence 
of a ferromagnetic-paramagnetic phase transition one usually restricts the study, within 
this class of models, to those where the random exchange couplings are exclusively fer- 
romagnetic. This restriction may seem simplistic, but in fact it turns out that many 
models exhibiting site or bond dilution (systems where some couplings are zero), or even 
asymmetric distributions of positive and negative couplings, belong to the same critical 
universality class. Sometimes these statistical models are called "weakly disordered" in 
order to distinguish them from the strongly disordered models encountered in spin-glass 
theory, but in fact this nomenclature is somehow inappropriate, since the disorder (whose 
strength is here given by the spread of the distribution of random- valued couplings) could 
either be weak or strong. It is often observed that models with varying disorder strength 
display identical critical properties, a situation reminiscent of the well-known universal- 
ity of the second-order phase transitions in non-disordered systems. Besides being of 
fundamental theoretical interest this universality is often of great practical importance: 
Whilst in calculations it is often simpler to consider weak randomness, in numerical stud- 
ies stonger randomness is preferred, in order to reach more easily the new critical regime 
induced by the disorder. 

An efficient laboratory for the study of disordered systems is provided by two-dimen- 
sional spin models. These models are non-trivial, critical, universal (unlike their one- 
dimensional counterparts) and easier to analyse, both analytically and numerically, than 
three-dimensional systems. 

In this broad playground many results have been obtained over the years. Analytic 
results were mainly derived from perturbative renormalisation group calculations [l], ||| 
whilst numerical results were obtained using Monte Carlo simulations and numerical 
diagonalisation of transfer matrices P, ^, . Much less has been found out concerning 
the conformal field theories (CFTs) which should describe exactly the random systems 
at their critical points. This is deceiving, considering the success of CFT in describing 
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the critical behaviour of the corresponding pure models (without disorder). Our study 
of coupled Potts models is mainly motivated by our interest to progress towards the 
identification of these conformal theories. 

The relation between the two problems, coupled and random, becomes clear when one 
considers the replica approach to random systems. For instance, the partition function 
of Ising-type models with random nearest-neighbour spin couplings can be cast in the 
form 

Z (P) = ex P \P X! JxiaVxVx+a \ = H eX P \ ~P H J x,a£x,a \ ■ (1-1) 
{a} I x,a J { CT } I x,a ) 

We can consider, for simplicity, the model on a square lattice. Then x runs over the 
sites of the lattice and the unit vector a over the neighbours, a = e 1; e 2 . The couplings 
J x>a assigned to links of the lattice take random values, independently for each link, with 
some specified distribution. The product of neighbouring spins —cr x a x+a in Eq. fll.l]) is 
the energy operator e XjCt assigned to a link, which becomes a local energy operator e(x) 
in the continuum limit near the critical point. In this limit the partition function fll.lf ) 
can symbolically be represented as 

Z(/3) oc Tr exp j-ft - J d 2 xm(x)e(x)^ , (1.2) 

Tio representing the Hamiltonian (or Euclidian action, in the field theoretic terminology) 
of the corresponding pure CFT if it is known. 

m ( x ) oc -r- , (1.3) 

is a mass-type coupling, which replaces (3J X)0t in the continuum limit: the randomness of 
J Xt a, whose average value is Jo, translates into the randomness of m(x). The trace Tr in 
Eq. ( |1.3| ) is assumed to represent, in the continuum limit theory, the summation over the 
spin configurations in the lattice model (|1 . 1| ) . 

In quenched disordered systems, averaging over the randomness should be done at 
the level of the free energy F(/3) oc log Z(f3). To this end one introduces replicas, that is, 
N copies of the same system 



( N N 

(Z(P)) N ck Tr exp - £ - j d 2 xm(x) £ e a { 

K a=l a=l 



x) (1.4) 



and averages over the randomness: 

WPW- (1-5) 



By taking the limit N — > one recovers the average of the free energy for a single system: 



[Z{(3)Y - 1 



F((3) cx log Z (13) cx hm . (1.6) 

If the distribution of the couplings { J x , a } in Eq. ( |1.1|) is gaussian, and consequently that 
of m(x) in Eq. ( |1.4| ), the disorder-averaged (Z(j3)) N in ( |1.5| ) is equivalent to a system of 
N models coupled by their energy operators: 

-£ft a) -m J d 2 xJ2^(x)+g I d 2 xY,£a(x)e b (x)\. (1.7) 

a=l a=l a^tb ) 

Here g$ is the width of the distribution of m(x), and m is its average value. At the 
critical point mo should vanish. 

In cases where the distribution of {J x , a }, or of m(x), is not gaussian, the resulting 
theory in Eq. (\1.7j ) contains higher-order coupling terms, like 

E e a (x)e b (x)e c (x) and ^ e a (x)e b (x)e c (x)e d (x), (1.8) 

a,b,c a,b,c,d 

but these are either irrelevant (in the renormalisation group sense) or do not modify 
the fixed point structure and its stability. We can thus limit our study to second order 
energy-energy couplings. 

On the lattice, before taking the continuum limit, a similar procedure applied to Z(f3) 
in Eq. (|1.1| ) produces 

(W= E^p(-^oEE4i + ^oEE4:^s|, (i-9) 

{<J a } [ X > a 1=1 X '' Q J 

where el a l = — cr| a Va+ a , and A is some constant coefficient. 



Thus, in order to solve the random coupling problem, one first has to solve the theory 
of iV coupled homogeneous models, Eqs. ( |1.7|) and ( |1.9|) , ultimately taking the limit 
N — > 0. A complete solution to the problem would be the identification of the exact 
conformal field theory associated with the critical point. 

When looking for possible candidates to this conformal theory, an important issue 
arising is the non-unitarity of the N — >• limit theory, as can be seen using perturbative 
renormalisation group (RG). For instance, in a non- unitary theory the central charge 
increases along the RG flow, in contradistinction to what is the case for a unitary (i. e. with 



reflection positivity) theory 0. As a first step, in order to avoid the problems of non- 
unitarity and work with a well-defined problem, we suggest to consider the problem of 
TV (positive, integer) coupled models, Eqs. (|1.7|) and ( |1.9|) , and to examine the critical 
properties of such systems. In the conformal theory language this amounts to studying 
N minimal models, coupled two by two by the operators $^2 which are the energy 
operators in the case of Ising or Potts models. This problem is unitary. The /5-function 
of the renormalisation group, to second order in perturbation, is given by 

% = (3{g) oc (N - 2)g 2 + 0(g 3 ) (1.10) 

in the case of Ising Jj]] , and 

^ = (3(g) oc 3eg + (N - 2)g 2 + 0{e\ g 3 ) (1.11) 

in the case of g-state Potts models @, ||, where e oc q — 2. By the usual analysis, for 
N < 2 (and eventually for N —>■ 0) the theory (|1.10 ) runs into zero coupling along the 
RG flow, and the theory (|1.11| ) runs into a new non-trivial fixed point with 

g FP = g* = ^ + 0(e 2 ). (1.12) 

Thus, the random Potts model possesses a new non-trivial fixed point, and it is therefore 
of interest to look for the associated conformal theory. 

These considerations hold true when g is initially positive, that is go > 0, which is 
the case for the random model. If now, in order to attain unitarity, we replace N < 2 
by iV > 2 (N integer), the coefficients in Eqs. fll.10 ) and (|1 . 1 1| ) change sign and the 
theory runs into a strong coupling regime, which is not controlled by pertubative RG 
and believed to be massive. This means a finite correlation length, and thus a non- 
critical theory. 

To avoid this problem, when passing from iV < 2 to N > 2 one should simultaneously 
change the sign of the initial coupling g Q . We expect some sort of similarity, or duality, 
of the domains N < 2, g > and N > 2, g < 0. In fact, if multiplied by (2 — N), 
Eq. Ql.ll ) takes the form: 

^ oc3eA- A 2 + C(A 3 ), (1.13) 



4 



where we have defined a new coupling A = {2 — N)g. In other words, the RG flows depend 
on g in the combination (2 — N)g, so that the regions iV < 2, g > and N > 2, g < 
are similar. 

In this way one can gain unitarity whilst still staying critical. Note that the critical 
coupling g*, Eq. ( |1.12p , depends only on N and e(q). It is not obvious, but we can hope 



that the unitary N > 2 problem being solved, that is the exact conformal field theory 
being found for general N, one could analytically continue into the domain N < 2 and 
extract exact information about the random model, N — ► 0. 

Following this approach, we first have to find the exact conformal theory of the non- 
trivial fixed point g* of N coupled Potts models, N — 3,4,5, .... The first step is to study 
three coupled models, expecting to be able to extend the result to general N later. It 
should be observed at this point that the model of N coupled Potts lattices (minimal 
conformal theories) is interesting on its own right, independently of its relation to the 
random problem. This makes the project doubly interesting. 

The system composed of three coupled state Potts models can be studied by the 
usual methods, in particular numerically, either by Monte Carlo simulations or by di- 
agonalising the transfer matrix in a strip geometry. Our first objective is to get some 
confidence in the existence of the fixed points predicted by perturbative RG, and second, 
to get fairly accurate numerical values for the central charge and dimensions of operators 
like spin, energy and their moments a\{x)a2{x), a\{x)a2{x)a^{x) , etc. This should be 
useful when searching for the corresponding conformal theory. 

To put the numerical study on a firm basis we need a properly defined model on 
the lattice, or rather on three lattices which are coupled, link to link, by their energy 
operators; see Eq. ( |1.9|) . For the Potts model the energy operator is of the form 

£x,a — ~&<j x ,<j x+a i (1-14) 

where S^y = 1 for o = a' and otherwise, and the spin operator a x at the lattice site x 
takes q different values. The partition function of three coupled models is of the form 



X ex P\f3 XX <W + + <W) + 9o XX'W^.t' + ^t,t'<W + <W<W) ] (1- 

{<x,t,??} V x,a x,a ) 



15) 



with a = a x and a' = a x+a etc. In this way we have simply copied Eq. (O) for the 



particular case of iV = 3 and a,b = 1, 2, 3, absorbing the coupling constant J into a 
redefinition of (3. 



Next we have to look for critical points of Eq. ( 1.15 ), for some (3 and negative 
If the position of the critical point is not known analytically and has to be determined 
numerically for a system as complicated as three coupled Potts lattices, with q 3 degrees 
of freedom at each site, accuracy is likely to be very poor, leading to unprecise or absurd 
measurements of critical properties. Thus, our numerical studies would greatly benefit 
from an exact determination of the position of the critical point. This can be done if the 
model is self-dual, like the single Ising or Potts models. For three coupled models, the 
existence of duality relations requires the inclusion of a three-energy interaction term 

x,a 

In the continuum limit this becomes 

d 2 x Si{x)s2{x)e?X x )- (1.17) 



It is straightforwardly shown that such a term does not modify the fixed point structure 
nor its stability. This means that adding ( 1.16Q to the lattice Hamiltonian should not 



modify the critical properies of the model, the continuum limit theory being the same. 
In this way we finally arrive at a lattice model with partition function 

H exp I l a (8a,* 1 + S T y + 8 v tf) + b (5 a y5 T y + 5 T y5 v ^ + 5^8^) + c5 a ^8 T y8 ViV >\ J 

{<T,T,r]} \X,Q! / 

(1.18) 

The coupling constants, a, b and c, can be chosen so as to render this model self-dual. 
In Sec. [| we shall establish the corresponding duality transformations. It will turn out 
that the model possesses not a point, but a line in parameter space (a, b, c) on which it is 
self-dual. On symmetry grounds we should expect that its critical points belong to this 
line. In fact, it will be shown that generically there exists three self-dual fixed points: 

1. That of a single Potts model with q 3 states of the spin variable. At this point 
a = b = and c ^ 0. Whenever the energy-energy coupling between the models is 
relevant (i.e., for q > 2), we have g 3 3> 4. The phase transition is thus first order, 
and the model is not critical. 
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2. That of three decoupled g-state Potts models. At this point b = c = and a ^ 0. 
The phase transition is second order if q < 4. 

3. The non-trivial fixed point of three coupled models. The transition is here second 
order, and the model is critical. The study of this point is the principal subject of 
this paper. 

To make a better connection with the continuum limit theory, the energy operators 
of the lattice model will be redefined in the next section to take the form 

E a)a r = l-6 aja ,. (1.19) 

The principal coupling term then becomes 

— 9o ^2 (Ev,a'E T ,T> + E TjT tE v rf + E^ a iE n rf) , (1.20) 

where the parameter go is a linear combination of the parameters b, c in Eq. ( [L.18 ): 

9o = b + c. (1.21) 

With respect to g the critical point |l|), corresponding to a first order transition, has 
positive g , as it should have been expected. The decoupling point |^), has g = 0, whilst 
the fixed point is found for finite negative go. 

As will be evident from the subsequent sections, to locate this last fixed point we 
have relied on the c-theorem of Zamolodchikov M, which states that the effective central 
charge of the theory decreases along the RG flow. The effective central charge has been 
measured using the strip geometry, as will be explained in detail in Sec. [|. We assumed 
on symmetry grounds, and verified numerically, that the RG flow from the decoupling 
point to the non-trivial coupling point goes along the line of self-duality that we have 
found. To locate the point |3|) we stayed on the self-dual line, on the negative go side, 
and followed, using transfer matrices on a strip, the evolution of the effective central 
charge along the line. This led us to a particular limiting point on the line of self-duality, 
actually its end-point, at which the exact values of the couplings are known. We then 
used transfer matrices and Monte Carlo simulations, with the couplings tuned exactly to 
their end-point values, to check scaling laws and measure critical dimensions of various 
operators, comparing the result with those obtained by perturbative CFT. 
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The paper is laid out as follows. In section 0, we present duality transformations and 
identify the self-dual lines for two and three coupled models. Section |3] is devoted to the 
computation of the central charge and the critical exponents of physical operators. To this 
end we employ perturbative CFT techniques. Section |] introduces the various transfer 
matrix algorithms we used to numerically compute the critical properties. Since these 
algorithms are interesting on their own right, some of them being more efficient than those 
previously described in the Potts model literature, they are presented extensively. Section 
[5] presents the numerical results obtained using these newly introduced algorithms. A 
Monte Carlo study of scaling laws is also undertaken. Finally, section ^| is devoted to 
concluding remarks and to a brief summary of the obtained results. 



2 Self-duality and criticality 

Duality relations, i.e. maps between sets of coupling constants that lead to the same 
partition functions, are central objects in the study of critical systems. These relations 
map one part of the coupling phase space to another one, a self-dual manifold separating 
the two. If the fixed points we are looking for are unique, they should be self-dual points 
in the phase space. If this were not the case, then fixed points would arise in dual pairs, 
which implies dual RG flows. Either these flows would never cross the self-dual line or 
they would both cross it at a given point. We are not able to picture in which way our 
system could behave like that. Numerical results presented later in this paper strongly 
support the conjecture that the fixed points of interest, for our models, are self-dual. 



2.1 Duality relations 

Let us first outline the construction of a global duality transformation from which the 
self-dual lines will be extracted. Following Ref. [|10J , we shall restrict our attention to 
Hamiltonians of the form 

n = -J2J($\---^P), (2.i) 

where = 1, . . . ,q are the spin variables of the kth model and Q,- = \Q k ' — 
Energy-energy coupled Potts models naturally have Hamiltonians of this form. For 
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generic J, the partition function can be written as 



Z = E II «(&). ( 2 - 2 ) 



where £jj = (Q 1 , • • -d^ ) an d the local Boltzmann weights are 

u$ = exp (/3J(D) • (2.3) 

The matrix [/, whose elements are w^(^), is a g x g cyclic matrix. It is rather obvious 
that the partition function can be rewritten as 

^9^114), (2-4) 

where a factor of has been pulled out in front, since after setting all the £y-s, one still 
needs to fix the absolute value of one spin in each model in order to completely specify 
the configuration. Using Fourier transform the partition function can be defined on the 
dual lattice. Since it is translationally invariant, the eigenvalues of the matrix U are 
given by 

m = E^(^^ju(£), (2.5) 
which implies, using inverse Fourier transform, that 

«(&') =J2 T (^V)HV)T*(^,V), (2.6) 

■n 

where 

T(Cv) = q^e W \^^\ . (2.7) 



This can be inserted in Eq. ( ^2|) , and leads directly to 

Z = ± Um,), (2-8) 

%J=1 (hj) 

where Nr> is the number of sites on the dual lattice. The additional factor of q comes 
from the fact that once all the fjij are set, one still has to define the absolute value on 
one site. On a self-dual lattice, such as the square, the number of spin sites on the direct 
lattice and on its dual are equal (neglecting boundary effects), and one thus gets 

Z(u) = qZ{\). (2.9) 

The transformations ( |2.5| ) are well-defined global duality transformations. From these, 



self-dual solutions can be extracted. Let us treat the cases of two and three coupled 
models separately. 
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2.1.1 Two models 



Consider the case of two coupled models, with Hamiltonian 

n = Y. H Hi ( 2 - 10 ) 

(id) 

Hij = -a (S ai ^. + 8 TitTj ) - b 5 ai!aj 5 Ti!Tr (2.11) 

The choice of sign for the coupling b conforms with existing computations for the random 
model. 

The duality relations, mapping the couplings (a, b) to (a*, b*), are given by 
e 2a* +fe * e 2a+fe + 2(g-l)e a + (g-l) 2 

e 2a+fe _ 2e a + 1 ' \ ■ ' 

e 2a+fe + (g - 2)e° - (g - 1) 

e 2a+6 _ 2e « + 1 ' 1 ' J 

The denominator is introduced to ensure that configurations of zero energy remain as 
such under duality. The self-duality condition constrains couplings a and b to satisfy 

e' = ^1^. (2-14) 

Two points of direct physical interpretation can be found on this line. First, the 
decoupled point 6 = for which a = ln(l + ^/g), which is the usual critical temperature 
P c (q) of the decoupled models. The second, at a = and b = ln(l + q), corresponds to 
a g 2 -state Potts model, for which f3 c (q 2 ) = ln(l + q). We remark that the case of two 



coupled models was previously considered by Domany and Riedel [11]. 



2.1.2 Three models 

The Hamiltonian associated with the case of three coupled models is 

H=J2n i j J (2.15) 

(id) 

- c S aiiaj 8 TitTj 8 m)Vj . (2.16) 
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The introduction of a three-coupling term is necessary to produce self-dual solutions, 
since it is generated by applying the duality relations to our original Hamiltonian. The 
duality transformations are given by 

3a * +3b » +c * e 3a+36 +C + 3{g _ 1)e 2a +b + 3^ _ + {g _ ^3 

1 e 3a+3b+c _ 3g2a+fe _|_ 3 g a _ J ' v ' > 

, +b , e 3a+36 +C + {2q _ 3)e 2a +b + (g 2 _ 4g + tya _ {g _ 1} 2 

e 3a+36+c _ 3 e 2a+6 _|_ 3 e a _ ]^ ' l Z ' iC 7 

e 3a+3b+c + ( g _ 3 ) e 2a+fe + (3 _ 2g ) g a + (g - 1) 
_ e 3a+36+c _ 3g2a+6 _|_ 3 g a _ \ ' v ' / 

Self-duality solutions are found to be 

b (2 + yg)e a -(Vg + l) 



c 



e 2a 



eC _ C 3 a 3(e a -l)(yg + l) + gyg + l 

(e*(2 + v /g)-(l + v /g))3 " ^ 

The two trivial points found in the two-models case also arise here; the decoupling point, 
for which b = c = and a = ln(l + y/q), and the g 3 -state Potts model, with a = b = 
and c = f3 c (q 3 ) = ln(l + g 3 / 2 ). 



2.2 A convenient reparametrisation 

There is a certain arbitrariness in the way we choose to parametrise a point on the self- 
dual line. Ideally, the couplings entering the lattice Hamiltonian would in some way 
be comparable to those of the field theory describing its continuum limit. Evidently, 
decoupling points should be associated with null couplings in both schemes, but this still 
leaves plenty of room for a reparametrisation. We propose such a reparametrisation which 
makes closer contact with the physical considerations put forward in the Introduction, 
and, at the same time, facilitates the implementation of numerical methods. 



2.2.1 General case 

Introducing the parameter x (x G [l,+oo[) in Eq. (|2.16| ), the self-duality relations can 
be rewritten as 

e" = ^ (2.2!) 
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eC = 3(^ + l)M* -T) 

(V5 + 2)* (i-l) 3 ^ ' 

where 7 = Since the decoupled models we study have ferromagnetic ground states, 
it is more convenient to trade the 5 au<7j for the operators 

l-r,.r. 1 4,,, (2.24) 

Doing so, the Hamiltonian becomes 

Uij = -(3a + 3b + c) + (a + 2b + c)(E ah>7j +E Ti , Tj + E ViiVj ) 
- (b + c) (E a . :a .E Ti:T . + E <Ju<T .E r]ur) . + E Ti:T .E VitV .) 
+ c E ai , aj E n , Tj E Vi:Vj . (2.25) 

The constant term, being the ground state energy, can be gauged out. Now define the 
new one- and two-energy coupling constants 

(3 = a + 2b + c, 

g = b + c. (2.26) 

With this set of couplings the Hamiltonian is turned into 

Uij = (3 (E aua . + E TiiT . + E vuv .) 

~ 9 {E (7i ,(T J E TitTj + E ai ^ ] E r]i ^ i] + E T ^ T] E m ^^) 

+ c E <T . :a .E TijT .E Tli:r) .. (2.27) 

On the self-dual line the couplings (5 and g are parametrised by 



-H 



jq + 2 x-1 



(2.28) 



_ 3(yg + l)x(x- 7 ) , , 

A similar reparametrisation for the case of two models leads to the Hamiltonian 

JUj = P(E ai:a . + E TitTj ) - g{E c ^.E TuT .), (2.30) 
12 



with 

2x 2 + x{q - 1)' 



e"' 3 



2x + (q-l) , 
e 9 = ^ 2.31 

ar 

This particular parametrisation has some definite advantages over the original one. 
First, all Boltzmann weights remain finite along the self-dual line, something useful both 
in numerical studies and in a comparison with perturbative CFT. Second, as we shall see 
in Sec. along the self-dual line, the coupling constant g has a sign in agreeement with 
perturbative computations. Finally, in the three-models case, the three-coupling term 
becomes infinite when x — > oo, leading to a null Boltzmann weight. This implies some 
simplifications of numerical studies, as we now show. 

2.2.2 Limits on the self-dual lines 

As we just said, the limit x — > oo will be of special interest in the following sections. There 
we can simplify numerical computations (both using transfer matrices and Monte Carlo 
simulations) by suppressing some of the Boltzmann weights attached to the links. The 
local Boltzmann weight W(Lij) for the degrees of freedom coupling spin sites i and j is 
completely determined by specifying the number Ly of layers having different spin values 
at the two ends of the bond (ij). The total Boltzmann weight of a spin configuration is 
then given by 

W = I] W(Lij). (2.32) 

<i,j> 

If we consider the case of two coupled models, Eq. (|2.30| ), three possible weights W(L) 
arise: 

W(0) = 1, W(l) = e-P, W{2) = e~ w+g . (2.33) 
Using the self-dual parametrised solutions, we see that in the limit x — > oo this becomes 

W(0) = 1, W(l) = 1/2, W{2) = 0. (2.34) 

For the case of three models, Eq. ( p.25|) , there are four possible weights 



W(0) = 1, W(l) = e-P, W{2) = e" 2/3+9 , W{3) = e" 3 ^ 3 ^. (2.35) 
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At x — > oo this becomes 

" f(0) = 1 - = sc^TJ- W ' (2) = 3(^TT? w '< 3 > = °' < 2 ' 36 > 

We shall use this limit extensively later on. The fact that configurations with one or 
more equal to three are forbidden greatly simplifies the problem. It even gives hope 
to solve exactly the lattice model at the critical point. In Sec. 4.5 below we shall see that 
it is possible to reformulate the Hamiltonians (p. 30 ) and ( j2.25| ) so as to obtain even more 



null Boltzmann weights, by passing on to a random cluster picture. 

3 Renormalisation group study of 
coupled Potts models 

Coupled Potts models have already been studied in details using renormalisation group 
techniques. Details can be found in references [0, |3], [L2|, [13], [T3|, |j"5| . We shall only present 
here a summary of these results, which are exposed extensively in the given references. 

The continuum limit of the models under consideration is defined by 

N „ N 



H = ^T / Ha >i - g J d 2 x^2ei(x)ej(x), (3.1) 



i=l ' ijtj 

where Si{x) is the continuum limit of E ai}(7 . [see Eq. (|2.24j) l and corresponds to the energy 
operator, and 7io,?. are the Hamiltonians of the decoupled Potts models in the continuum 
limit. We shall restrict the discussion to g-state Potts models with 2 < q < 4, for which 
the dimension of the energy operator varies between 1 (for q = 2) and 1/2 (for q — 4). 
For such values of q, the continuum limit of the models are conformal field theories. The 
term g J ee is considered as a perturbation, which is relevant by dimensional analysis 
(except for q = 2, where it is marginal). The computation of the /3-function to two loops 
was done in Refs. @, and is given by 

P(g) = = 3eg(r) + 4vr(iV - 2)g\r) - 16vr 2 (iV - 2)g\r) + 0{g\r)). (3.2) 

a log(r) 

Here e is the perturbation parameter, and it corresponds to the dimension of the pertur- 
bation: e = | (1 — A £ ) where A e is the dimension of the energy operator of the decoupled 
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models. For q = 2 we have e = 0, whilst for q = 3, e = 2/15 and for q = 4, e = 1/3. The 
parameter r is an infrared cut-off. 

In order to compute critical exponents we need to compute the matrices {Z £ )ij and 
(Z a )ij which are respectively the renormalisation constants of the energy and spin oper- 
ators. These have to be understood in the following way: 

6 R = (Z )0, (3.3) 

where the ith. element of vector O is the operator O of model i, and Or is the correspond- 
ing renormalised operator. These matrices are, again to second order in the perturbation, 
given by @ 

rfl ° g( f; (r))tJ = (-MN - 1)9 - 8vr 2 (iV - I),? 2 + 0(g 3 )) (1 - <%) , (3.4) 
dl0g{Z d ° {r)) v = (3(N - l)gVeF + A(N - 1)(N - 2)n 3 g 3 + 0{g*)) (3.5) 

with 

T ~ 2 !*(-£)!*(-!) ■ (3 - 6) 



The fact that the energy matrix is not diagonal was observed in Refs. 0, [15| . This implies 



that the energy operators for each individual layer are no longer the proper ones to study 
critical behaviour. Instead, the eigenvectors of the matrix turn out to be the ones we 
observe in numerical studies. We shall come back to this point when we compute critical 
exponents in the next sub-section. 



The possible behaviours of our coupled systems can be summarised as follows |T3 |. 



For q = 2 the model corresponds to the iV-colour Ashkin- Teller model. For iV = 
2 it is integrable [T(|. For N > 2 the sign of g is determinant for the large- 
scale behaviour: When g > a second order phase transition of the Ising type is 
observed, whilst for g < the scenario is that of a fluctuation-driven first order 
phase transition |T7|, M. 



For q > 2 and N = 2 the model is still integrable, but now presents a mass 



generation [19|. This again indicates a first order phase transition. 
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• For the case N > 2 and g > the coupling constants flow far from our perturbative 
region. Even if a definite proof is not given, a comparison with the case q = 2 seems 
to tell us that a mass gap is dynamically generated, indicating a first order phase 
transition. 

The situation is completely different for g < 0. In that case there is a non-trivial 
infrared fixed point for 

The critical exponents associated with the energy and spin operators for this fixed 
point, along with the values of the central charge, are computed below. 

3.1 Critical exponents for q > 2, N > 2 

We now compute the critical exponents of our coupled models at the fixed points of the 
renormalisation group. For the decoupling point, the critical exponents are evidently the 
ones of the pure models. We shall therefore concentrate on the non-trivial fixed point 
identified above. Since the renormalisation matrix for the spin operator is diagonal and 
the coupling is invariant under permutation, any linear combination of the spin operators 
has a well-defined critical exponent which is found to be || 

where T was defined in Eq. ( |3.6|) . 

Defining critical exponents for the energy operators is somehow more tricky. Correla- 
tion functions between renormalised operators (at a given cut-off R) are of the following 
form 

(e,(0) £j -(i2)> ~ EE( z e)*Wii^<^(0)ei(l)>, (3.9) 

mixing correlation functions of the N different layers. One way to extract unambiguous 
critical exponents is to diagonalise the renormalisation matrix, thus introducing a new 
basis of operators. This diagonalisation is exact, in the sense that eigenvectors for the 
one-loop renormalisation matrix remain eigenvectors to all order in the perturbation. For 
three models the eigenvectors are 

e 1 + e 2 + e 3 , ei-e 2 , e 2 -e 3 . (3.10) 
16 



When using this basis, the computation of the different exponents is straightforward. We 
have JTH 

A £l+£2+£3 = A £ + 6e - 9e 2 + 0(e 3 ) (3.11) 

and 



A £l _ £2 = A £ -3e + -e 2 + 0(e 3 ). (3.12) 



9 
— ( 

2 

Some explicit values of the different critical exponents for the case of three coupled 
model are provided by Table |l|. The appearance of negative magnetic exponents for 
sufficiently large q shows the limitations of the perturbative expansion. It is also possible 
to compute the critical exponents for higher moments of the spin and energy operators. 
The physically significant operators are found to be 

cricr 2 , (J1CX3, <7 2 cr 3 , 0-1O-2CT3; (3.13) 
eie 2 + e 2 e 3 + e 3 e 1} eie 2 - e 2 e 3 , e 1 e 2 -e 3 ei, e x e 2 e z . (3.14) 

As was shown for the random Potts model [pOfl , perturbative computations for higher 
moments of the spin and energy operators are much less precise than for the operators 
themselves, and one should keep in mind that they eventually become absurd for suffi- 
ciently high moments (for example, the third magnetic moment has negative exponent 
for q > 3.7). For the spin operators, the second order perturbative computations lead to 
the following exponents PT], ^U] 

A CTlff2 = 2A CT (e) + (l - ((N - 2) log2 + ^)) + 0(e 3 ), (3.15) 



4(JV-2) V iV-2 V v J 12 
A aia2a3 = 3A CT ( e ) + — ^— (l - ( (N - 2) log2 + H + " ) ) + Q(e% (3.16) 



aia2C3 ayi A(N-2)\ N-2V ' 6 12 24 
with 

a = 33 — . (3.17) 

For the energy operators, critical exponents for the diagonal operators are given by 
(eg = e x e 2 + e 2 e 3 + e 3 e x and e\ = e 1 e 2 - e 2 e 3 ) 

A £| = 2A £ (e) + 3e+^ 2 + 0(e 3 ) (3.18) 

A 4 = 2A £ (e) - \e - 9e 2 + 0(e 3 ) (3.19) 

A £l£2£3 = 3A £ (e)-| e 2 + 0(e 3 ) (3.20) 
The numerical values for these operators are given in Table 2. 
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Q 






A CT1 


2.00 


1.0000 


1.0000 


0.12500 


2.25 


1.1447 


0.8499 


0.12789 


2.50 


1.2639 


0.7154 


0.12964 


2.75 


1.3615 


0.5930 


0.12985 


3.00 


1.4400 


0.4800 


0.12805 


3.25 


1.5006 


0.3737 


0.12353 


3.50 


1.5429 


0.2710 


0.11501 


3.75 


1.5624 


0.1656 


0.09926 


4.00 


1.5000 


0.0000 


0.04238 



Table 1: Perturbative critical exponents of physically significant energy and spin opera- 
tors for three coupled g-state Potts models. 



Q 










Acri CT20-3 


2.00 


2.000 


2.000 


3.000 


0.2500 


0.3750 


2.25 


2.005 


1.834 


2.837 


0.2775 


0.4553 


2.50 


2.021 


1.653 


2.664 


0.2949 


0.5190 


2.75 


2.046 


1.456 


2.479 


0.3030 


0.5685 


3.00 


2.080 


1.240 


2.280 


0.3023 


0.6048 


3.25 


2.126 


0.997 


2.060 


0.2921 


0.6283 


3.50 


2.186 


0.713 


1.806 


0.2703 


0.6375 


3.75 


2.272 


0.350 


1.486 


0.2303 


0.6268 


4.00 


2.500 


-0.500 


0.750 


0.0975 


0.5301 



Table 2: Perturbative critical exponents of physically significant energy and spin operator 
moments for three coupled g-state Potts models. 
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3.2 Central charge and the c-theorem for q > 2, N > 2 



One of the most important tools of CFT, Zamolodchikov's c-theorem ||, provides a 
simple way of computing the central charge for perturbed theories. Moreover, it gives us 
crucial information: 

• The c-function, to be defined below, is decreasing along the renormalisation flow. 

• If the flow has a fixed point, the field theory at the fixed point is conformally 
invariant and its central charge is the value of the c-function at that point. 

With our conventions for the /3-function taken into account, the c-function is defined 



as 



c(g') = c pmc -2Aj o g /3(g)dg, (3.21) 

with Cpure being the total central charge of the decoupled models. The central charge 
deviation from the decoupling point value is thus easily computed. Substituting g' = g*, 
the fixed point value of the couplings, we get the following correction: 

Ac = -24 P /3(g)dg (3.22) 



_ 27 N(N-l) f 9 A 5 

~ 8 (N-2)* { 2{N-2) € ) +U{€) [6 - Z6) 

Table ^ presents the pure and perturbative values of the central charge. Clearly, for large 

enough q the latter are not to be trusted, since the perturbation theory eventually breaks 

down. This is witnessed by the change of sign in Eq. (|3.23|) when e > |. 



4 The transfer matrices 

Systems of several coupled g-state Potts models possess an enormous number of states, 
and one should think that accurate numerical results cannot be obtained by the transfer 
matrix technique, since only very narrow strips can be accessed. However, there are 
several possibilities for dramatically reducing the size of the state space, as we show in 
the present Section. In particular we shall see that it is possible to numerically study the 
g-state Potts model for any real q with less computational effort than is required in the 
traditional transfer matrix approach to the Ising model. 

19 



Q 


Cpure 


CFP 


2.00 


1.5000 


1.5000 


2.25 


1.7627 


1.7620 


2.50 


1.9975 


1.9931 


2.75 


2.2089 


2.1976 


3.00 


2.4000 


2.3808 


3.25 


2.5734 


2.5500 


3.50 


2.7309 


2.7164 


3.75 


2.8734 


2.9054 


4.00 


3.0000 


3.3750 



Table 3: Comparison of the central charge values at the pure and the non-trivial fixed 
points. 

We have optimised the algorithm described in Ref. ||, and adapted it to the case of 
three coupled models, by taking into account all known symmetries of the Boltzmann 
weights in the spin basis. Further progress can be made by trading the spin variables for 



clusters or loops. The cluster algorithm of Ref. |22[ has been adapted to the problem of 
coupled models, and we also describe, for the first time in the literature, the practical 
implementation of the associated loop algorithm. The latter algorithm allows us to 



address the case of three coupled models on strips of width L r 



12. 



4.1 The four algorithms 

Consider a system of coupled Potts models defined on a cylinder of circumference L and 
length M, measured in units of the lattice constant. The imposition of periodic boundary 
conditions in the transversal direction is understood throughout. The transfer matrix can 
be viewed as a linear operator T that acts on the partition function of the M row system, 
where fixed boundary conditions have been imposed on the degrees of freedom of the Mth 
row, so as to generate the corresponding quantity for a system with M + 1 rows. To fix 
terminology, we shall refer to the specification of the boundary condition on the last row 
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as the state of the system. Thus, symbolically, 

where the Greek subscripts label the possible states. The matrix elements T aj g are nothing 
but the Boltzmann factors exp[— 7i(a, (3)] arising from the interaction between the degrees 
of freedom in the Mth and the (M + l)th row, when these are in the fixed states (3 and 
a respectively. 

It is well-known that crucial information about the free energy, the central charge, and 
the scaling dimensions of various physical operators can be extracted from the asymptotic 
scaling (with strip width L) of the eigenvalue spectrum of the transfer matrix, and so 
the question arises what is the most efficient way of diagonalising it. A quite common 
trick is to decompose it as a product of sparse matrices which each add a single degree 
of freedom to the Mth row, rather than an entire new row at once. But even more 
important is the realisation that the total number of states, and thus the dimensionality 
of the transfer matrix, can be reduced by taking into account the various symmetries of 
the interactions and possibly by rewriting the model in terms of new degrees of freedom 
other than simply the Potts spins. We shall refer to any such collection of states as a 
basis of the transfer matrix, and the cardinality of the basis as its size. The accuracy 
of the finite-size data, of course, improves significantly with increasing L, and since the 
size is a very rapidly (typically exponentially) increasing function of L we must face the 
highly non-trivial algorithmic problem of identifying and implementing the most efficient 
basis. 

These considerations have led to the development of four different algorithms for the 
coupled Potts models problem. We list them here in order of increasing efficiency and, 
roughly, increasing algorithmic ingenuity. Details on the implementations will be given in 
the following sub-sections. As a rough measure of their respective efficiency we indicate 
how large L can be attained in practiced] for the three-layered model. 

• algl uses the trivial basis of Potts spins and can access L = 5 for q = 3 and 

L = 4 for q = 4. Although hopelessly inefficient in the general case this algorithm 
3 The largest sparse matrices that we could numerically diagonalise using a reasonable amount of 
computation time had size ~ 10 7 . 
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is still of some use for small, integer q since it allows for a straightforward access 
to magnetic properties. 

• alg2 still works in the spin basis, but the Z q symmetry of the Potts spins has been 
factored out. In this process the magnetic properties are lost. The size is now 
independent of q, but for q < L a further reduction of the size occurs due to a 
truncation of the admissible states. For q = 3 this algorithm can access L — 7, but 
for general q only L = 6. 

• alg3 utilises a mapping to the random cluster model, and the basis consists of the 
topologically allowed ('well-nested') connectivities with respect to clusters of spins 
that are in the same Potts state. Magnetic properties are lost, but can be restored 
through the inclusion of a ghost site at the expense of increasing the number of 
basis states. Since q enters only as a (continuous) parameter, this algorithm is par- 
ticularly convenient for comparing with analytic results obtained by perturbatively 
expanding in powers of (q — 2). In the non- magnetic sector L = 7 can be accessed. 

• alg4 works in a mixed representation of random clusters and their surrounding 
loops on the medial lattice. Magnetic properties can be addressed through the 
imposition of twisted boundary conditions, which by duality are equivalent to the 
topological constraint that clusters and loops do not wrap around the cylinder. 
Again q is a continuous parameter. Due to the definition of the medial lattice, only 
even L are allowed. This algorithm can access L = 12 in the non-magnetic sector 
and L = 10 in the magnetic one. 

Before turning our attention to a detailed description of these algorithms we briefly 
recall how physically interesting quantities may be extracted from the transfer matrix 
spectra. 

4.2 Extraction of physical quantities 

The reduced free energy per spin in the limit M — > oo of an infinitely long cylinder is 
given by 
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where Ao is the largest eigenvalue of T. Starting from some arbitrary inititial vector of 
unit norm |v ), this can be found by simply iterating the transfer matrix 



Ao(L)= lim i-l„||T^|vo)||. (4.3) 

Higher eigenvalues A^(L) are found by iterating a set of n vectors {|v fc )}^~Q, where a 
given |vfc) is orthogonalised to the set {(v;)}^ 1 after each multiplication by T |24 . 



In general, of course, there exists more expedient methods for diagonalising a square 
matrix, but since each of our four algorithms allows for factorising the transfer matrix as 
a product of sparse matrices this simple iteration method is superior. The advantage of 
using sparse-matrix factorisation is that it reduces the number of elementary operations 
from (Cl) 2 to LCl per iterated row pjfl , where Cl is the size of the basis. 



It is well-known from conformal field theory how to relate the central charge to the 



finite-size scaling of the free energy [26 



TTC 

/oW=/o(oo) - — + •■•. (4.4) 

Similarly, the gaps of the eigenvalue spectrum fix the scaling dimensions Xi of physical 
operators^ [^8j 

fi (L)-f (L) = ^ + .... (4.5) 

In general one may construct several sectors of the transfer matrix by identifying 
the various irreducible representations of the symmetry group of the microscopic (spin) 
degrees of freedom. As a familiar example consider just one single Ising model (q = 2). 
Here there are two sectors (even and odd) corresponding to the possible transformations 
of a state vector |v) — > ±|v) under a global spin flip 

a — > (a + 1) mod q . (4.6) 

This generalises straightforwardly to the g-state Potts model for which we must consider 
transformations |v) — > e 2i7r 9 |v), j = 0, 1, . . . , q — 1. The thermal and magnetic scaling 



dimensions, x t and then extracted by applying Eq. (|4.5|) to f° vcn — /q v611 and 

/o dd — /o vcn respectively. In the ordinary spin basis (as used by algl) the choice between 
4 The notational discrepancy with Sec. |^ is intentional. Indeed, the identification between the scaling 
dimensions of physical operators and numerically measured quantities is the object of Sec. ||. 
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the even and odd sector can be very easily accomplished by iterating initial vectors 
(cfr. Eq. pop ) that are either even or odd under the transformation ( |4.6| ). When the 
Z q symmetry has been factored out (as in alg3 and alg4) the odd sector has to be 
constructed explicitly by means of a ghost site or a seam (see below for details). 

When considering several coupled Potts models more than one magnetic exponent 
may be defined. Namely, for each individual model one may independently choose be- 
tween the even (e) and odd (o) sector. At the decoupling point, of course, one has 

but at a general critical fixed point the magnetic exponents thus defined are independent. 

When using Eqs. ( |4.4| ) and ( [4.5| ) to obtain finite-size estimates for c and Xi the conver- 
gence properties can be considerably improved by explicitly including higher-order terms 
in 1/L. In an extensive study of the g-state Potts model Blote and Nightingale [29] nu- 



merically showed that the sub-leading correction to the free energy takes the form of an 
additional 1/L 4 term on the right-hand side of Eq. ( |4.4|) . Not surprisingly an analogous 
statement can be made about Eq. ( |4.5| ) for the scaling dimensions. From the viewpoint 
of conformal field theory such a non-universal correction to scaling can be rationalised as 
arising from the operator TT, where T denotes the stress-energy tensor |5| . This operator 
of dimension 4 is present in any theory that is conformally invariant pC]. 

Finite-size estimates c(L, V) can then be extracted either from two-point fits for 
fo(L) and fo(L' = L + 1) using Eq. ( [4.4|) or from three-point fits for fo(L), f (L + 1) and 
f (L' = L + 2) using 

/o(i) = /oN-^ + | + -, (4.8) 



Similarly, for the scaling dimensions the one-point estimates Xi(L) obtained from Eq. (4~5) 
can be improved by extracting two-point estimates Xi(L, L + 1) from fits of the form 

/i(i)-/o(i)~ + | + - (4-9) 



24 





g = 2 




q = 3 




q = 4 




N 


2 


3 


2 


3 


2 


3 


L = 2 


16 


64 


81 


729 


256 


4,096 


L = 3 


64 


512 


729 


19,683 


4,096 


262,144 


L = 4 


256 


4,096 


6,561 


531,441 


65,536 


16,777,216 


L = 5 


1,024 


32,768 


59,049 


14,348,907 


1,048,576 




L = 6 


4,096 


262,144 


531,441 




16,777,216 




L = 7 


16,384 


2,097,152 


4,782,969 








L = 8 


65,536 


16,777,216 










L = 9 


262,144 













Table 4: Number ^ (L) of basis states in algl as a function of strip width L and the 
number N of coupled Potts models. We indicate here these numbers up to the largest 
strip width for which we were able to diagonalise the transfer matrix. 

4.3 Algorithm algl 

In the trivial spin basis each state is specified by the values of the L Potts spins in row 
M. The size of the transfer matrix for N coupled models is therefore 

Kn( L ) = 9 LN - (4-10) 

These numbers increase very rapidly as a function of strip width L, and they do not 
have a g-independent upper bound. To highlight the merits of the more sophisticated 
algorithms that we are about to develop we present some explicit values of the N^ 1 in 
Table |. 

The reason that the trivial spin basis is still of some use is that it does not explicitly 
break the Z q symmetry of the Potts spins. Therefore, T contains both the even and 
the odd sectors, and the various magnetic scaling exponents can easily be extracted. 
Furthermore, algl has served as a check of the results obtained from the optimised 
algorithms presented below. 

In the case of the Ising model (q = 2), states which are even and odd under a global 
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spin flip [cfr. Eq. (16)] are given by 



| e > = -L(|00---0>+|ll •••!>), |o) = -^(|00-.-0>-|ll..-l>). (4.11) 

This is easily generalised to several layers of spins by simply forming direct product 
states. For example, for two coupled models we can define the states 

|eo) = |e) ® |o), |oo) = |o) ® |o), (4.12) 

so that the gap between the largest eigenvalue in each of these sectors and the largest 
eigenvalue in the totally symmetric (|ee) = |e) ® |e)) sector defines two magnetic expo- 
nents, x e £ and which are in general independent. 

For the g-state Potts model the appropriate prescription reads 
|e> = i-(|00...0> + |ll-..l) + ...+ |( g -l)( g -l)...( g -l))), (4.13) 

|o) = -L(\00...0) + e 2 T |ll...l) + ... + e 2 ^|(g-l)(g-l)...(g-l)) 
where the physical state is obtained from Eq. (|4.12|) by taking the real part. 



4.4 Algorithm alg2 

When the Z q symmetry of the Potts spins is factored out the size of the basis can be 
dramatically reduced. The drawback of this approach, however, is that the magnetic 
properties are lost. In other words, the resulting transfer matrix does only have an even 
sector, but this is still enough to extract finite-size estimates for the central charge and 
the thermal exponent. (On the other hand, we shall see in the following sub-sections 
that in the case of the algorithms alg3 and alg4 it is possible explicitly to reconstruct 
the odd sector from topological considerations.) 

The algorithm alg2 was already used in Ref. []5J to study the random-bond Potts 
model, and before adapting it to the case of several coupled Potts models we briefly 
recall its application to the single- layered model. 

The basic idea is that the 5^,0- -type interactions do not depend on the explicit values 
of the spins ai and Oj. What matters is whether the two spins take identical or different 
values. Therefore, the possible number of states for a row of L spins is equal to the 
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number of ways bi in which L objects can be partitioned into indistinguishable parts 
| ST| . With m v parts of v objects each {y = 1, 2, . . .) this can be rewritten as 

oo oo r I 

= E ' n T^n- (4-14) 

m„=0 t/=l 

where the primed summation is constrained by the condition Y^=\ vm^ = L. Alterna- 
tively one can write || 

2 m3 m,4 

6l=EEE-E1. (4-15) 

«2=1«3 = 1«4 = 1 «L=1 

with ink = 1 + max(i2, «3, • • • , ifc-i). From the former representation the generating 
function is found as 

exp(e< - 1) = £ (4.16) 



whence 



n=0 n! 



1 00 £* L 



Yet another way of interpreting the b^ is to notice || that they are the total number 
of possible L-point connectivities, including the non-well nested ones [29] . We shall come 
back to the notion of well-nestedness when we discuss alg3. 

A major advantage of the basis used in alg2 as compared to the trivial spin basis is 
that the bi do not depend on q. Thus, any integer value of q can be accessed with this 
algorithm. However, for any fixed q the size of the transfer matrices that have L < q can 
be further reduced.0 Namely, in this case the number of permissible states is truncated 
due to the fact that L objects cannot be partitioned into more than L indistinguishable 
parts! In this way we identify the number of states N^ 2 (L) for N = 1 Potts layers as a 



sum over Stirling numbers of the second kind [[32 



Nff{L) = £ st ] = e A E(-ir~ fe TV • ( 4 - 18 ) 



m=l m=l m ' k=0 



Explicit values for q = 2, 3, 4 are shown in Table and we conjecture that asymptotically 
N i 2 (L) ~ q L . We recall that the upper limit for general q is bi, and this increases as 
L aL , where a is a constant of order unity 0. However, by comparing Tables [| and |5] we 
see that even though the number of states in algorithms algl and alg2 exhibit the same 
asymptotic growth, alg2 is much more efficient than algl. 
5 This possibility was not discussed in Ref. ||. 
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9 = 2 






q = 3 






q = 4 






N 


1 


2 


3 


1 


2 


3 


1 


2 


3 


L = 2 


2 


4 


4 


2 


3 


4 


2 


3 


4 


L = 3 


4 


10 


20 


5 


15 


35 


5 


15 


35 


L = 4 


8 


36 


120 


14 


105 


560 


15 


120 


680 


L = 5 


16 


136 


816 


41 


861 


12,341 


51 


1,326 


23,426 


L = 6 


32 


528 


5,984 


122 


7,503 


310,124 


187 


17,578 


1,107,414 


L = 7 


64 


2,080 


45,760 


365 


66,795 


8,171,255 


715 


255,970 




L = 8 


128 


8,256 


357,760 


1,094 


598,965 




2,795 


3,907,410 




L = 9 


256 


32,896 


2,829,056 















Table 5: Number N^ 2 (L) of basis states in alg2. For iV = 2,3 coupled models these 
numbers are shown up to the largest strip width L for which we were able to diagonalise 
the transfer matrix. 

4.4.1 Layer indistinguishability 

In the general case of N coupled Potts models further progress can be made by observing 
that since all layers interact symmetrically there is an additional permutational sym- 
metry. If we imagine numbering the one- layer states by an integer i = 1, 2, ... , N q ^(L) 
a general iV-layer state can be represented by (h,i2, ■ ■ ■ , ijv) where i\ > > . . . > ijy. 
The total number of states for iV coupled models is then 

= EE £•••£!■ (4-i9) 

11=1 12=1 13=1 *JV = 1 

For iV = 2, 3 this is easily found to be 

N q ,s = i (JVJ X + 3iVj x + 2iV ftl ) . (4.20) 
Again explicit values pertaining to alg2 are shown in Table |5|. 
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4.5 Mapping to random cluster model 



An altogether different approach for setting up the Potts model transfer matrices consists 
in rewriting the partition function in terms of extended objects (clusters and loops) rather 
than the local spins. The resulting random cluster models and loop gases |J3J have 
the advantage that the specification of their states no longer depends on the value of 
q. Instead q enters only through the fugacity of the non-local objects, and hence it can 
be taken as a continuously varying parameter. Such reformulations of the problem are 
especially convenient for making contact with the predictions of conformal field theory, 
and in particular with the perturbative expansion in powers of (q — 2) which has been 
studied in Sec. |[ 

The practical implementation of transfer matrices for such non-local degrees of free- 



dom was pioneered by Blote and collaborators [^2, |35|. The basic idea is here that in 
two dimensions the allowed connectivities of clusters and loops are strongly constrained 
by topological considerations. Accordingly the size of the corresponding transfer matrix 
is drastically reduced. 

In the following two sub-sections we generalise these algorithms to the case of iV 
coupled Potts models, and we give explicit details on the implementation for iV = 1,2,3. 
The cluster algorithm alg3 has previously been described for N = 1 by Blote and 
Nightingale |[29|| . It was discussed in more detail in Ref. ||, where it was also shown 
how to adapt it to the case of bond randomness. The loop algorithm alg4 has to our 
knowledge not previously been used to study the Potts model.f] It is however the most 
efficient algorithm that we know of, and in fact it performs much faster than the spin 
basis algorithm for the Ising model! 

Before focusing on the concrete implementation of alg3 and alg4, which is the object 
of the following two sub-sections, we dedicate this sub-section to developing the appro- 
priate mappings of the partition function. We need to consider the cases of two and 
three coupled Potts models in turn, but it should be clear that the results generalise 
straightforwardly to any number of models N. 



'Details on the implementation of related but more complicated loop models can be found in Refs. BE 
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4.5.1 Two models 



The partition function for two coupled Potts models can be written as 

z = e n ex p(^)> ( 4 - 21 ) 

W^} (v) 

where Hij = a(<5 <TiiCT . + 8 TitT .) + b5 au(T .8 Ti!T .. Using the standard Fortuin-Kasteleyn trick 
p3| the exponential of Hij is turned into 

exp(Hij) = (1 + u a S ff .^){l + u a 5 n>Tj )(l + u b S aii(7j 5 Ti!T:j ), (4.22) 

where we have defined u a = e a — 1 and tt& = e fe — 1. After some straightforward algebra 
we obtain 

exp(Hij) =k + k 1 (5 u ^ (Jj + 5 n:Tj ) + k 2 8 au(Tj 8 TuTp (4.23) 

where 

k = l, h = u a , k 2 = u 2 a + u b (l + u a ) 2 . (4.24) 



Now imagine expanding the YIua product in Eq. ( |4.21|) , using Eq. ( [4.231) . In this 
way we obtain a total of 3 E terms (E = \{ij) \ being the number of lattice edges), each 
consisting of a product of E factors of k{ [i — 0, 1, 2) multiplying a product of Kronecker 
deltas. Define £ to be the graph consisting of two copies ('layers') of the lattice on which 
the Potts model is defined, one copy placed on top of the other. To each of the terms 
in the expansion we associate a graphical representation Q on £ by colouring the lattice 
edges for which the corresponding Kronecker deltas occur in the product. In other words, 
Q takes the form of a two-layered bond percolation graph. 

To reproduce the partition function (|4.21|) we need to perform the sum J2{a,r} over the 



spin degrees of freedom. Because of the Kronecker deltas this is trivially done for each 
term Q, yielding simply a factor of q for each of the C connected components ('clusters') 
in (?.[] The factors of fcj are easily accounted for by writing E = E + E\ + E 2 , where Ej 
is the number of occurences in Q of a situation where precisely j = 0,1,2 edges placed 
on top of one another have been simultaneously coloured. The partition function then 
takes the form 

Z = Y j q C kf'k^kf. (4.25) 
Q 



7 Note that a single isolated site is to be counted as a cluster on its own. 
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On the self-dual line ( 2.14Q the edge weights ( [4.24 ) assume the simple form 



k = l, h = e a -l, k 2 = q. (4.26) 

As has already been mentioned, the limit a —>■ oo is of special interest. In this limit 
we can rewrite the partition function as Z = exp(Ea)Z , where Z has the same form as 
( [4.25| ) but with the finite edge weights 

k = 0, k\ = 1, k 2 = 0. (4.27) 

It is worthwhile noticing that in this limit it suffices to specify the colouring configuration 
of one of the layers to deduce that of the other layer. Namely, whenever an edge in the 
first layer is coloured its counterpart in the second layer has to be uncoloured, and vice 
versa. Therefore, the configuration sum J2g appearing in Eq. ( |4.25| ) in fact runs over one 
layer only. 

4.5.2 Three models 

In the case of three coupled Potts models the nearest-neighbour interactions take the 
form TLij = a5i + b5 2 + c5 3 , where we have introduced the short-hand notation 

Si = S ai!<Tj + 5 TuTj + 8 vum , 

h = 6 aii(Tj 5 Ti7Tj + 5^6^ + K,Tj$m,vv ( 4 - 28 ) 

S3 = ^ i ,<j j 8 T .,r j K,m- ( 4 - 29 ) 
The objects Si, 5 2 , S3 are easily shown to obey the following relations 



<5i<5i — Si + 25 2 , 5i5 2 = 25 2 + 35 3 , 6x63 = 3S 3 , 
S 2 S 2 = 5 2 + 6S 3 , 5 2 5 3 = 35 3 , 5 3 5 3 = 5 3 . 



(4.30) 

Simple algebraic manipulations then lead to the following result, generalising Eq. ( 4.23| ). 



ex.p(Hij) = k + k ± 5i + k 2 5 2 + k 3 5 3 (4.31) 

with edge weights 
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u„ 



k 2 = u 2 a + u b (l + u a f , 

h = ul + 2>u a u h {l + u a f + (1 + u a f \u c (l + u b ) 3 + 3u 2 b + u\ 



(4.32) 



We recall that the Fortuin-Kasteleyn parameters are w a ,&, c = e a ' 6,c — 1. After the original 
spin degrees of freedom have been summed over, the partition function takes on the form 



g 



(4.33) 



with a notation analogous to that employed for the case of two coupled models. The 
graph configurations Q now consist of three bond percolation graphs stacked on top of 
one another. 



Along the self-dual line (|2.20| ) the edge weights ( |4.32| ) simplify to 

ft = l, h = e a -l, k 2 = qV 2 (e a -l), h = q 3 ' 2 . 



(4.34) 



Finally, in the limit a — > oo the modified partition function Z = exp(—Ea)Z has the 
finite edge weights 



k = 0, k 



1 = a, 



k 2 = q l/ \ h = 0. 



(4.35) 



4.6 Algorithm alg3 

With the mapping to a random cluster model in hand we are now ready to discuss the 
implementation of alg3. The point of depart of this algorithm is the form ( f4.33| ) [or 



( |4.25| )1 of the partition function, in which the spin degrees of freedom have been turned 
into non-local clusters. 

For simplicity we consider first the case of a single Potts model, with partition function 
Z = Y^g l ^ 1 [|33|| . To construct the transfer matrix we need to specify a basis, so that 
the knowledge of the state before and after the addition of a new degree of freedom gives 



us enough information to compute the appropriate Boltzmann weights, cfr. Eq. (f4.1|). As 



shown by Blote and Nightingale [^3] this is achieved by specifying the connectivity (with 
respect to the clusters in Q) of the L points in the last row of the strip. Since connections 
can only be mediated by the lattice edges which have previously been added there is 
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a very powerful topological constraint (known as well-nestedness) on the connectivity 
states. Namely, given any four consecutive points, if the first is connected to the third 
and the second to the fourth, then all four points must be connected. Consequently, the 
number cl of allowed (well-nested) L-point connectivities is less than the total number 
of connectivities bi considered in Section O. 



Using a recursive principle the well- nested connectivities can be enumerated |22|, and 
they turn out to be nothing but the Catalan numbers 

c L = zy Ml iyy ~ 4* for L » 1. (4.36) 

Since the number of Potts states enters only as a (continuous) parameter in the parti- 
tion function this is independent of q. A mapping from the connectivities to the set of 
consecutive integers 1, 2, . . . , cl can also be established. More details on the practical 
implementation of the transfer matrices can be found in Ref. ||. 

For N > 1 layers of coupled Potts models we need to simultaneously keep track of the 
connectivity state of each layer in order to compute the factors of q occuring in Eq. (|4.25 ) 



or ([4.33|) . To specify the state of the layered system we can however take advantage of the 
indistinguishability of the layers, cfr. Eq. ( |4.20| ). The resulting number of connectivity 
states N^ g3 (L) used in alg3 can be found in Table |^. 

4.6.1 Magnetic properties 

At the expense of increasing the size of the basis it is possible to generalise alg3 to 
treat the case of a Potts model in a uniform magnetic field H ||22|| . To this end, note 
that the interaction with the field can be accounted for through the inclusion of a term 
HJ2i^cri,cr in the Hamiltonian, where a so-called ghost spin of fixed value a® = 1 has 
been introduced. Taking each site of the lattice to be connected to through a 'ghost 
edge', this has the usual form of a nearest-neighbour interaction. The mapping to the 



random cluster model therefore goes through in exact analogy with Sec. [4.5| . In specifying 
the connectivities one now needs both to keep track of which sites are connected (directly 
or indirectly) to the ghost site and, at the same time, to specify how the remaining sites 
are interconnected. 
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L alg3 


£alg4 


iV = 1 


N = 2 


N = 3 


2 


2 


2 


3 


4 


3 


4 


5 


15 


35 


4 


6 


14 


105 


560 


5 


8 


42 


903 


13,244 


6 


10 


132 


8,778 


392,084 


7 


12 


429 


92,235 


13,251,095 


8 


14 


1,430 


1,023,165 




9 


16 


4,862 


11,821,953 





Table 6: Number iV^ ls3 ' 4 (L) of basis states in the even sector of algorithm alg3 and alg4. 
For N = 2, 3 coupled models these numbers are shown up to the largest strip width L for 
which we were able to diagonalise the transfer matrix. With alg4 it is possible to access 
larger strip widths L alg4 = 2(L alg3 — 1) than with alg3, using the same number of states. 
Since for a strip of width L, alg4 needs to employ intermediary states of (L + 2)-point 
connectivities, the size of the basis given here pertains to these intermediary states. 
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The extended L-point connectivity states thus defined can be ordered and enumerated 
\T1 1, and their number is found to grow asymptotically as 5 L . For H — > the transfer 
matrix has the following block form 



(4.37) 

where superscript 1 (0) refers to the (non-)ghost connectivities, and the magnetic scaling 
dimension xh is obtained via Eq. ( f4.5|) from the largest eigenvalues of the T 00 and T 11 
blocks. The physical content of this relation is that by acting repeatedly with T 11 on 
some initial (row) state |vo) 7^ one measures the decay of clusters extending back to 
row 0. This must have the same spatial dependence as the spin-spin correlation function 
and hence be related to xh f22fl . 

For L — 1, 2, 3, . . . the size of the magnetic block T 11 is d^—c^ = 1,3,10,37,146,599,. . . . 
For three coupled models this means that computations with one, two or all three layers 
in the magnetic sector are feasible for strip widths up to L = 6. Since our most refined 
algorithm alg4 can access L = 10 we turn our attention to this next. 




4.7 Algorithm alg4 

The configurations of the random cluster model on some graph are in a one-to-one corre- 
spondence with configurations of a fully-packed loop model on the medial graph |53J. By 
definition the nodes of the medial graph are situated at the mid-points of the edges of 
the original graph. In our case the graph on which the Potts model is defined is simply 
the square lattice, and the medial graph is then nothing but another square lattice that 
has been rotated through ir/4 with respect to the original one, and rescaled by a fac- 
tor of 1/ \[2. Bipartitioning the medial graph into even and odd sublattices the precise 
correspondence is as shown in Fig. |l]. 

The partition function for N coupled Potts models can now be rewritten in the loop 
picture by using Euler's relation p4^ . Namely, for each layer the number of clusters (C), 
loops (I), coloured edges (e) and sites (s) are related by 2C = 1 + s — e. Thus, in the case 
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even site 



odd site 



Figure 1: The relation between the random cluster model on the square lattice and 
the loop model on the corresponding medial graph. The clusters consist of connected 
components of coloured edges (thick lines) or isolated sites (filled circles). Loops on the 
medial graph (dashed lines) are defined by the convention that they wrap around the 
cluster boundaries. 



of two models, Eq. (|4.25 ) is turned into 

Ei / , \ E 2 



*=«-XW(;k) '(t)'- (4 ' 38) 



Q 

Note that the factor of q~ e has been redistributed using e = E\ + 2E 2 . Similarly, for the 
three-layered model Eq. ( |4.33| ) is transformed into 



g \y 



h \ E1 ( h \ E2 ( h ^ Ez 



(4.39) 



When constructing the transfer matrix the advantage of this representation is that 
less information is needed to keep track of the loops than was the case in the cluster 
representation. Roughly speaking this is because the loop model is fully packed, so that 
one does not need to waste information to specify where is 'the empty space' in between 
the clusters. 

The strip width L is now defined as the number of 'dangling ends' resulting from 
cutting the loops in between two rows of sites on the medial graph. A sparse-matrix 
decomposition can then be made by adding one vertex at a time, rather than one edge 
as was the case in the cluster representation. This has been illustrated in Fig. ||] for a 
strip of width L = 6. For each added vertex there are two possible configurations of the 
loop segments that must be summed over, cfr. Fig. [I]. When adding the first vertex of a 
new row the number of dangling ends increases from L to L + 2, and it only goes back 
to L once the full row of L vertices has been completed. Therefore, we need to work 
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interchangingly with bases specifying the loop connectivities of L and L + 2 dangling 
ends [SB] . 

By definition of the loops, each of the dangling ends is connected to exactly one other 
dangling end. In particular the strip width L must be even. Employing basically the 



same recursive argument as for the clusters P21 it is found that the possible number of 



connectivities among L dangling ends is now only ||35|| . For L ^> 1 this increases as 
2 L , and in fact the loop representation of the Potts model is even more efficient than the 
spin representation of the Ising modelQ This is witnessed by Table BJ, where we show 
some explicit values of N^ S \L) = N^ g3 ((L + 2)/2). 

At this point a brief comment on the boundary conditions is in order. Since the 
medial lattice is rotated through 7r/4 with respect to the original one, the imposition 
of periodic boundary conditions on the loops is not a priori equivalent to the boundary 
conditions hitherto used for the clusters. Indeed, these two possibilities are connected 
by a modular transformation, as should be clear from Fig. ^. The consistency between 
the results obtained from alg3 and alg4 will thus serve as a useful check of the modular 
invariance of the critical system under investigation. 

To implement Eqs. fl4.38| ) and ( |4.39[ ) for several coupled Potts models we need to keep 
track of the edge weights on the original lattice as well as the number of closed loops on 
the medial graph. As shown in Fig. [I] the loop configuration suffices to determine the 
positions of the coloured edges on the original lattice, and so the edge weights can easily 
be determined from each vertex appended to the medial graph. By the same token the 
single-layer algorithm furnishes an even more efficient way of performing transfer matrix 
calculations for the random-bond Potts model than the one presented in Ref. H. 



4.7.1 Magnetic properties 

In alg4 the magnetic sector of the transfer matrix is constructed by using the fact 
that the spin-spin correlator maps onto a disorder operator under duality 0. For a 

single Potts model at the self-dual (critical) point, or for several coupled Potts models 
8 Using Stirling's formula a more accurate estimation of the asymptotic behaviour of c L j 2 is found 
as 2 L (L + 2) _3 / 2 ^/8/7r. This approximation is asymptotically exact in a strict sense, and its relative 
precision is better than 10 % for L > 6. 
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along the self-dual lines described in Sec. |2|, this leads to a relation between the magnetic 
scaling exponent and the largest eigenvalue of a modified transfer matrix on which twisted 
boundary conditions have been imposed. 

Defining the local order parameter as 

M a (r) = (^ a{rU - ^ , a=l,...,q. (4.40) 

it is easily seen that the magnetic two-point correlator G aa (ri,r 2 ) = (M a (ri)M a (r 2 )) is 
proportional to the probability that the points r\ and r 2 belong to the same cluster. Let 
us briefly recall that any given configuration of the random clusters is dual to one where 
each coloured edge in the direct lattice is intersected by an uncoloured edge in the dual 
lattice, and vice versa. Taking the two points ri and r 2 to reside at opposite ends of 
the cylinder, the graphs contributing to G aa (ri, r 2 ) thus correspond to dual graphs where 
clusters are forbidden to wrap around the cylinder. This is equivalent to computing the 
dual partition function with twisted boundary conditions 

a — > (a + 1) mod q (4.41) 

across a seam running from r\ to r 2 . By permuting the Potts spin states the shape of 
this seam can be deformed at will as long as it connects r\ and r 2 . 

A realisation of the Potts model transfer matrix in the presence of a seam was first 
described in Ref. || within the context of alg3. Since there is a one-to-one correspon- 
dence between clusters and loops that wrap around the cylinder these ideas can also be 
applied to alg4. 

In Fig. [| we show a configuration contributing to the partition function with twisted 
boundary conditions. As the square lattice is self-dual we shall take the clusters and loops 
to live on the direct lattice and its medial lattice respectively. The extended connectivity 
state of the L (resp. L + 2) dangling ends is the direct product of the cl/2 possible 
connectivities mediated by the loops, and an integer specifying the position of the seam. 
The seam is a path inside the infinite dual cluster spanning the length of the cylinder, 
and after the addition of each vertex it must be updated according to the invariant that 
no loop closure take place across the seam. 

Now consider the situation shown on Fig. |2| where we are about to add the fifth vertex 
of the top row. Since this is an even site the two possibilities for the loop configuration 
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Figure 2: Computation of the partition function with twisted boundary conditions. 
Coloured edges (solid linestyle) connect the sites of the direct lattice (filled circles) so as 
to form clusters. The loops on the medial lattice (short boldface dashes) surround the 
clusters, and both cluster and loops are subject to periodic boundary conditions (indi- 
cated by open circles) across the strip of width L = 6. A seam (long dashes) connecting 
sites of the dual lattice prevents the clusters and loops from wrapping around the cylin- 
der. By adding a single vertex at a time the transfer matrix can be decomposed as a 
product of sparse matrices. Each multiplication by a single-vertex matrix corresponds 
graphically to augmenting the 'jigsaw puzzle' by one of the two 'pieces' shown in the top. 
In the situation at hand, one of these is forbidden by the twisted boundary conditions. 
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L 


A/mag 1 


N — 2 

1 v mag ^ 


N — 3 

J ' mag 


2 


12 


20 


20 


4 


225 


600 


680 


6 


5,880 


22,344 


30,856 


8 


189,630 


930,510 


1,565,620 


10 


6,952,176 


41,451,696 


83,112,744 



Table 7: Number (L) of basis states needed by alg4 to access the scaling exponents 

of the three-layered system that correspond to the insertion of a magnetisation operator 
in iV ma g = 1, 2, 3 of the layers. 

are those shown in the left part of Fig. |l[ However, the first of these would lead to a 
forbidden loop closure, as witnessed by the fact that the seam is 'trapped' and cannot 
be updated. The corresponding entry in the transfer matrix is therefore zero, and only 
the graph realising the second possibility contributes. 

Using the numbering of the L (resp. L + 2) dangling ends of an (in)complete row 
shown at the bottom (top) part of Fig. ^ it is easily seen that the seam position is always 
immediately to the right of an even end in every other row, and to the right of an odd end 
in the remaining rows. Therefore, in all cases the number of permissible seam positions 
is equal to half the number of dangling ends. The total number of extended L-point 
connectivity states is therefore 

Lc L/2 _ L\ 



2 n . 



(4.42) 



2 (f -!)!(* + 1)1 

Taking into account that for a strip of width L we need intermediate states of (L + 2)- 
point connectivities, the size A^ s of the transfer matrix for three coupled models with 



mag 



1,2,3 layers in the odd sector is as shown in Table U\. Note that for N, 



mag 



= 1,2 
3 all 



two of the layers are indistinguishable in the sense of Sec. |4.4.1| , whilst for iV mag 
three layers are indistinguishable. 

A very important remark pertains to the proper choice of the initial vector |v ) used 
for finding the largest eigenvalue, cfr. Eq. ( |4.3|) . Namely, when more than one layer is in 
the odd sector the seam positions on all layers must coincide at points r\ and r 2 . The 
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reason why this is so can readily be illustrated for the case of two coupled models. With 
the initial vector (symbolically) chosen as |v ) = |1, 1, 1, . . .) <g> |1, 1, 1, . . .) one would 
observe the asymptotic scaling of the correlator 

( E M£\r 1 )M£\r\)M£\r 2 )Mf){r> 2 )\ ~ ( £ M a (n)M a (r 2 )\ , (4.43) 

\r 1 ,r' 1 ,r 2 ,r' 2 I \n,r2 I 

where the summations run over the L points at each extremity of the cylinder. At large 
distances |t"2 — r*i] ^> 1 one would expect this to scale with the exponent 2x^\ where 
Xft is just the scaling dimension of the usual magnetisation operator. However, with 
|v ) = |1, 0, 0, . . .) <8> |1, 0, 0, . . .) one would instead observe the scaling of 

(E MW(n)MW ( ri )MW(r 2 )MW(r 2 )) , (AAA) 

and this should decay as |r 2 — ri\~ 2x( n\ where x$ is the sought scaling dimension of 
the local operator M^\r)M^(r). Indeed this is the simplest example of a multiscaling 
exponent as discussed by Ludwig |2j in the context of the random-bond Potts model 

The proper choice of |v ) is tantamount to anchoring all the seams at the point 
r\. When applying Eq. fl4.3|) one should then theoretically project Tl r2 ~ ri l|v ) out on 
a state with a definite and identical seam position for all the layers before taking the 
norm. However, since we expect all the non-zero entries of the iterated vector to grow as 
exp (|r 2 — ri|A(L)) this would just correspond to multiplying by a constant before taking 
the logarithm. Clearly such a constant would not contribute to the computed value of 
A(L), and so the projection can be omitted. 



5 Numerical results 

Using the transfer matrix algorithms just described we are able to find the effective values 
of the central charge and of the various thermal and magnetic scaling dimensions, all as 
functions of the coupling constant a, parametrising the position on the self-dual lines 
identified in Sec. |2|. We are furthermore able to monitor the scale dependence of these 
quantities by changing the strip width L. 
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Let us briefly recapitulate the physical information that we hope to extract from 
these data. First, we should like to provide compelling evidence that, for N > 2 coupled 
Potts models, a novel, non-trivial critical fixed point is located on the self-dual line, in 
the limit a — > oo. We shall presently describe how such a conclusion may be attained 
from our numerical data. Second, we aim at fixing the values of the critical exponents 
at this fixed point. This is done by taking the limit a — > oo explicitly in the transfer 
matrices, in order to obtain high-precision data for rather large strip widths (up to 
L — 12 for N — 3). Extrapolating these data to the infinite system limit L —>■ oo we 
shall be able to identify the a — > oo fixed point with the one obtained from perturbative 
CFT in Sec. |3|, and to associate the numerically obtained scaling dimensions with those 
of physical operators in the continuum limit. Third, we provide an independent check 
of the criticality of the fixed point under consideration. This is done by verifying the 
existence of scaling laws, using Monte Carlo simulations directly in the limit a — > oo. 
As a by-product we extract values of the scaling dimensions that agree with the more 
precise estimates obtained from the transfer matrices. And fourth, we use our transfer 
matrices to inquire further into the structure of the (presently unknown) CFT governing 
the critical behaviour of the a — > oo model. To this end we take a closer look at the 
higher eigenvalues in the even sector. These data determine the scaling dimensions of 
less relevant operators in the Verma module, and they give crucial information of the 
operator content and descendence structure of the sought CFT. Finally, a comparison 
of the results obtained from algorithms alg3 and alg4 serves as a check of the modular 
invariance of the theory. 

Although our main interest is N = 3 coupled models we have also produced numerical 
results for the case of two coupled g-state Potts models, for several values of q G [2,4]. 
For q = 2 this is the Ashkin- Teller model. This model is critical on the entire half line 
a G [0, oo], and it provides a useful check of our algorithms since the critical exponents 



are known exactly as a function of a [37], |34|. For q > 2, on the other hand, Vaysburd 
|19| has predicted a dynamical mass generation and thus a non-critical model. Since this 
prediction was made under several assumptions it is reassuring to see that our numerics 
confirms it. At the same time, our N = 2 data provide a clear illustration of the difference 
between a first order phase transition scenario and a second order one. 
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5.1 Finding critical points from effective exponents 



The key to the search for critical fixed points is, of course, Zamolodchikov's c-theorem 
||. This theorem states that the effective central charge (the c-function) decreases along 
the RG flow and is stationary at its fixed points. In view of our definition ( p. 21 ) of the 



c-function, this is equivalent to the familiar statement that the /3-function vanishes at a 
fixed point. 

A local minimum (maximum) of the effective central charge thus corresponds to a 
(un)stable fixed point. When, as is the case here, the RG flow takes place in a multi- 
dimensional space of coupling constants one can also imagine the existence of saddle 
points, corresponding to a partially stable fixed point. But generically we expect the 
flow to take place along the self-dual line, and accordingly we can limit the search to 
this one-dimensional self-dual submanifold. This assumption is nicely corraborated by 
our numerical results. Indeed we find that in general any motion perpendicular to the 



self-dual line leads to a decrease of the central charge [see, e.g., Fig. y.b and Fig. [10 
below]. Hence, invoking duality, this line must serve as a mountain ridge for the central 
charge. 

A notable exception to this scenario happens for the Ashkin- Teller model, where the 
self-dual half line a G [0, oo] bifurcates into two mutually dual lines of critical points at 
a = 0. However, this anomaly is very clearly signalled by the vanishing of the c-function 
of the other self-dual half line, a G [—oo, 0[. We have not observed a similar situation for 
any other value of q, or of N. 

Once a candidate for a fixed point has been identified the question arises whether it 
is critical or not. To address this point it is useful to examine the dependence of c e g on 
the system size L. If a critical point is involved, the finite-size estimates c e s(L) should 
converge very rapidly, with increasing L, to the true value of the central charge at this 
point. On the other hand, if the point is not critical a finite correlation length £ must 



be involved. As has been shown by Blote and Nightingale |29| for the particular case of 
a first order phase transition, Eq. Q4.4| ) for the finite-size scaling of the free energy must 
then replaced by 

e -£/£ 

fo(L) = /o(oo) - const.— — H . (5.1) 
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This in turn implies that the effective central charge, i.e., the one measured assuming the 
scaling form (|4.4j ), will vanish exponentially as a function of the strip width L. We shall 
soon see that our numerical data distinguish very easily between these two situations. 

Complementary information about the location of possible fixed points is provided 
by the effective values of the thermal and magnetic scaling dimensions x e s(L) along the 
self-dual line. In many situations a fixed point is signalled by the fact that the curves 
x e ff(L) and x e s(L + 1) intersect at some coupling a(L). The finite-size estimates a(L) 
then converge very rapidly towards the true value of the fixed-point coupling, as L — ► oo. 
Indeed, this is nothing but the well-known technique of phenomenological renormalisation 
p8| . Evidently a fixed point can also be characterised as a point of high symmetry. It is 
thus hardly surprising to see from our numerics, that in many cases local extrema of x e g 
as a function of a serve equally well to locate the position of the fixed point. 

As was the case for the central charge, once a fixed point has been located its precise 
nature (critical or not) can be inferred by observing the size- dependence of x e g (L). Since 
a finite correlation length implies an exponential, rather than a power law, decay of the 
various correlators, effective values of the scaling dimensions extracted from Eq. (|4.5| ) 
will be observed to vanish exponentially with L, if £ < oo. 



5.2 Two coupled models 
5.2.1 Ashkin- Teller model 

The case of two coupled Ising models is nothing but the well-known Ashkin- Teller model 
PP| . We shall nevertheless begin by investigating this case in order to demonstrate the 
method outlined above, and to see how we can reproduce known results. 



As has been shown by Baxter 33] the isotropic Ashkin- Teller model can be mapped 
to a staggered eight- vertex model. This fact that this model is staggered seems to impede 
its solvability, but on the self-dual line given by Eq. (|2.14|) , a £ [—oo, oo], the eight- vertex 



model becomes homogeneous and hence soluble. For a > the model is critical, with the 
following exact values of the critical exponents J3"7|: 



1, x t- 2 - y > ^"S' X "- 8 _ % , (5-2) 
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Figure 3: Central charge for two coupled Ising models, (a) shows the effective central 
charge as a function of the parameter a. (b) corresponds to moving perpendicularly off 
the self-dual line. 



where 



2/1.1 

y = -arccos ( e~ a + -e~ 2a 

y 7T V 2 2 



(5.3) 



In Fig. [3|.a we show the numerical values for the effective central charge, in the form of 



the three-point fits, cfr. Eq. (|4.8|) . For every a < the central charge tends to zero in the 
large-system limit, indicating non-critical behaviour, and for a > it approaches unity, 
in complete agreement with Eq. ( |5.2| ). The transfer matrix has also been diagonalised 
directly at a = oo, again with the result c ~ 1. 

Fig. |3|.b depicts the behaviour of the central charge along a line perpendicular to the 
self-dual line at the point a = 2. The value a = 2 is chosen arbitrary, and we observe 
similar curves for other values of a > 0. As we move into the non-self dual regime the 
central charge decreases, and the larger the strip width the faster the decrease. In the 
L — » oo limit we should have c c g = exactly. We can therefore conclude that the self- 
dual line indeed acts as a ridge of the central charge, thus confirming the hypothesis that 
the RG flow is along that line. We also observed that for a < the critical line bifurcates 
into two mutually dual Ising-like lines |34 . 



Next, in Fig. El, we present measured values of the thermal scaling dimension, x t , 
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Figure 4: Thermal dimension Xt and magnetic dimensions x^j , for two coupled Ising 
models. 

as well as the two magnetic ones, Xjj and x% . As the /3-function is exactly zero we 
can expect very accurate values, and as the comparison with the analytical results (|5.2| ) 
shows this is indeed the case. The only exception is the point a = 0, which corresponds 
to the 4-state Potts model. Due to the presence of a marginal operator there are strong 
logarithmic corrections |40|| , and the analytical results x t = 1/2 and x^ = x^ = 1/8 are 



not accurately reproduced. Finally, by taking the explicit a —>■ oo limit in the transfer 
matrices we have verified the results Xt = 3/2, x$ = 1/8 and xjj = 3/8. 

5.2.2 Two coupled 3-state Potts models 

In the case N = 2, q = 3 there is an a <-> —a symmetry at the level of the free energy, as 
is reflected by the values of the effective central charge shown in Fig. |^. The origin of this 
symmetry becomes clear if one considers the possible local Boltzmann weigths W(Lij, a) 
[cfr. Eq. ( |2.32| )1 associated with the interactions along edge (ij). For two coupled g-state 



Potts models these are given by 

W{Q,a) = l, W(l,a)=e a , W(2, a) = e 2a+b = 2e a + q - 1. (5.4) 

In the transfer matrix, starting from an arbitrary state, the first weight will appear 
(q — l) 2 times, the second one 2(q — 1) times, and the last one just once. Specialising 
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Figure 5: Central charge for two coupled 3-state Potts models (a) and two coupled 4-state 
Potts models (b). The dots in (b) correspond to an interpolaton from the measurements 
made directly at a = — oo and a = —5. We also indicate by a cross a fit at a = — oo for 
strip widths L=10-14. 

now to q = 3, we see that 1^(0, a) and W(l, a) have the same degeneracy, and that 

W(0,-a) = e- a W(l,a) } W(l, -a) = e~ a W(0, a), W(2, -a) = e~ a W{2, a). (5.5) 

Thus we have established an exact symmetry for the specific free energy on the self dual 
line 

/(-a,g = 3) = /(a,g = 3)-2o. (5.6) 

Returning to the central charge values shown in Fig. |awe see a local maximum at 
a = log(l + v3), the fixed point corresponding to two decoupled 3-state Potts models. 
At this point the finite-size estimates converge very rapidly towards c = 2 x |. The other 
trivial fixed point is located at a = and corresponds to a single 9-state Potts model. 
Here the estimates c e &{L) decrease steadily with increasing strip width, as expected 
for a first order phase transition. Any model with a bare coupling ao > log(l + \/3) 
renormalises towards the a — > oo fixed point, but this again corresponds to a first order 
transition [|l^] as witnessed by the steady decrease of the estimates for c. 
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L 


c(L,L + A) 


x t (L,L + 2) 


4 


2.251 


1.032 


6 


2.231 


0.856 


8 


2.203 


0.808 


10 


2.175 


0.762 


12 




0.720 



Table 8: Effective exponents for two coupled 4-state Potts models at a = — oo. 

5.2.3 Two coupled 4-state Potts models 

The effective central charge for two 4-state models can be inferred from Fig. |5].b, and 
apart from the disappearence of the a —a symmetry conclusions are as in the case of 
two 3-state models. The range of a-values shown on the figure cannot entirely exclude 
the possibility of a novel fixed point at a = — oo, but the data of Table |£| obtained by 
using alg4 directly at this point, tell us that if a = — oo indeed is a fixed point it must 
be unstable and non-critical. 

5.2.4 The point a = oo 

As has already been stressed in Sec. ^ the topological algorithms alg3 and alg4 enable 
us to treat q, the number of Potts spin states, as a continuously varying parameter. We 
finish the discussion of the two-models case by presenting some results for the effective 
exponents for several fractional values of q, obtained directly in the a — > oo limit. 

First, in Table |9| we display the central charge values. For q = 2 these converge very 
nicely towards c = 1, the exact Ashkin- Teller value ( |5.2| ). For higher q the increasing 
spacing between succesive estimates signals non-critical behaviour, or in other words, the 
free energies are not well fitted by Eq. ( f4.4|) . It is convenient to have such data at one's 
disposal, since they facilitate the distinguishing between first and second order phase 
transitions in the sequel. Evidently, for q = 2.25 it not a priori easy to clearly point out 
the first order nature of the transition, since the correlation length is very large. However, 
it is interesting to notice that whilst for q = 2 the estimates increase monotonically as a 
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q 


c(4,8) 


c(6, 10) 


c(8, 12) 


c(10,14) 


c(12,16) 


2.00 


0.988504 


0.995838 


0.998256 


0.999142 


0.999528 


2.25 


1.142687 


1.149660 


1.151260 


1.151271 


1.150843 


2.50 


1.263404 


1.264818 


1.260899 


1.255641 


1.250201 


2.75 


1.351006 


1.338262 


1.319944 


1.300188 


1.280109 


3.00 


1.406234 


1.367834 


1.322338 


1.274092 


1.224246 


3.25 


1.430583 


1.353605 


1.266530 


1.174515 


1.079463 


3.50 


1.426433 


1.298639 


1.157785 


1.011642 


0.864787 


3.75 


1.396976 


1.208954 


1.008126 


0.807868 




4.00 


1.346005 


1.092727 


0.833815 


0.590962 





Table 9: Effective central charge for two coupled g-state Potts models at a — > oo, obtained 
from alg4. 

function of L, for all q > 2 they eventually begin to decrease for large enough L. 

Table [10] provides the analogous estimates for the thermal scaling dimension. For 
q = 2 they converge towards x t = 3/2, in accordance with Eq. ( |5.2| ). Otherwise their 
variation is not monotonic in L, presumably signalling the presence of a finite correlation 
length. 

5.3 Three coupled models 

We now turn to our primary interest, namely the case of N = 3 coupled Potts models. 
The self-dual line here starts at the point where b = — oo, that is a e [a min , oo] with 

a min = log J . (5.7) 

Unlike the case of two models just discussed we shall find that the a — > oo limit cor- 
responds to a critical fixed point for all values of q G [2,4]. The determination of the 
associated critical exponents thus becomes of paramount interest, and we have dedicated 
considerable computational effort to this purpose. Let us discuss the various cases in 
turn. 
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Q 


zt(6,8) 


^(8,10) 


x*(10,12) 


z t (12,14) 


z t (14,16) 


2.00 


1.467238 


1.475693 


1.490527 


1.495490 


1.497566 


2.25 


1.745413 


1.554317 


1.588008 


1.605994 


1.618971 


2.50 


1.863461 


1.768649 


1.702778 


1.739814 


1.769780 


2.75 


1.780405 


1.717780 


1.672429 


1.635165 


1.603300 


3.00 


1.707436 


1.630297 


1.575072 


1.532502 


1.499100 


3.25 


1.648514 


1.565197 


1.510808 


1.474971 


1.452970 


3.50 


1.605054 


1.524930 


1.482309 


1.464823 


1.466231 


3.75 


1.576886 


1.508815 


1.487488 


1.497813 




4.00 


1.562949 


1.514505 


1.521838 


1.566459 





Table 10: Effective values of the thermal scaling dimension for two coupled g-state Potts 
models at a — > oo, obtained from alg4. 

5.3.1 Three coupled Ising models 

This is the three-colour Ashkin- Teller model. The effective central charge is displayed 
in Fig. |||.a. At the trivial fixed point at a = we find the usual dip in c, reflecting the 
first-order nature of the transition in the 8-state Potts model. 

For all a > the central charge is approximately constant, c ~ 3/2, and we have also 
verified this value in the limit a — > oo. This situation is very reminiscent of that of the 
(standard) two-colour Ashkin- Teller model, and seems to allow for a situation where the 
entire half line a > consists of critical fixed points along which the critical exponents 
vary continuously. Whether this is indeed the case will be the subject of a separate 
publication |41| . 



Fig. |6].b presents a closer look at the antiferro magnetic region a < 0. We find here 
a surprising candidate for an unstable critical fixed point at a ~ —0.22 with a central 
charge of c ~ 1.5. Further support for this suspicion is found in the maximum for the 
effective thermal exponent on Fig. |7].b, and in the crossing of the curves for the magnetic 
exponents on Fig. |].b. More results along these lines will be published elsewhere [j4l"H . 

The results for the thermal exponent on the half line a > are shown on Fig. [7|.a. 
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ure 6: Central charge for three coupled Ising models. 
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Figure 8: Magnetic exponents for three coupled Ising models. 



Everywhere the convergence seems to be very rapid, except at a = where we expect 
an extrapolated effective exponent of = due to the first order transition. In the 
a — > oo limit we estimate x t = 1.24 ± 0.01; more detailed results are provided by Table 
1~2] below. 

Finally, in Fig. |8| we present results for the three different magnetic exponents, x$ , 
Xff and Xff . The first exponent is constant, Xjj ~ 1/8 for a > 0, as in the case of two 
coupled Ising models, whilst the other two depend on a. In the a — > oo limit we find 
x$ ~ 0.27 and x^ ~ 0.46 as witnessed by Tables [T5| and |TJ. 



5.3.2 Three coupled 3-state Potts models 



Until now, the models for which we have presented numerical data are believed to be 
somewhat atypical representatives for the general class of iV coupled (/-state Potts models. 
On one hand, for N = 2 coupled models the 0(g 2 )-term in the /3-function Ql.ll ) vanishes, 
and neighter perturbative RG nor numerics predict any non-trivial critical fixed points 
for q 7^ 2. On the other, for q = 2 the energy-energy coupling is exactly marginal, leading 
either to the Ashkin- Teller scenario with an entire line of critical fixed points along which 
the exponents vary continuously, or simply to a flow back towards the decoupled fixed 
point. For the case of iV = 3 coupled 3-state Potts models, which we consider now, the 
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Figure 9: Central charge for three coupled 3-state Potts models. The dots designate an 
interpolation from a = 4.5 to a = +oo. 

situation can thus be believed to be generic: Perturbation theory predicts a novel critical 
fixed point, Eq. J1.12Q , with unique critical exponents, and so does numerics, as we shall 
presently see. 

Consider first the effective central charge along the self-dual line a G [a m i n , oo], where 
a min = log • Ori Fig. |9|.a we recognise the familiar structure of a local minimum 

and a local maximum, signaling the usual two trivial fixed points. However, for larger 
values of a a novel feature emerges, as witnessed by Fig. |9|.b. The central charge is here a 
very slightly decreasing function of a, indicating a flow towards the fixed point at a = 00. 
Furthermore, the finite-size dependence of the estimates leads us to the conclusion that 
this point is now critical. 

The more accurate data of Table [lT] corroborate this scenario. First, we see that the 
finite-size estimates converge very rapidly to the value 



c = 2.377 ±0.003, 



(5.8) 



in very good agreement with the perturbative prediction cfp = 2.3808 + 0(e 5 ) from 
Eq. (|3.23| ).P| Second, the convergence of the estimates in Table [11] is now from below, 



9 In fact, the latter prediction would be slightly smaller if we could go to even higher order in pertur- 
bation theory, since the series ( 3.23| ) is known to be alternating. 
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Figure 10: Central charge for three coupled 3-state Potts models when moving perpen- 
dicularly off the self-dual line. 

convincingly excluding the possibility that c e g could eventually tend to zero, as would 
be the case for a non-critical fixed point. Third, the numerical value for c is comfortably 
away from that of the decoupled fixed point, c pure = 3 x | = 2.4. 



In Fig. [10] we show the behaviour of the central charge along a line perpendicular 
to the self-dual line at the point a = 1.19217. Once again, this value of a was chosen 
arbitrarily; similar curves are observed for other values of a > 0. 

Plots of the effective thermal and magnetic exponents are shown in Figs. |IT].a and 



12. a respectively. More accurate values obtained directly at the a = oo critical point are 



given in Table |T^-|T6|. 

5.3.3 Three coupled 4-state Potts models 

The situation for three coupled Potts models is very similar to the case of q = 3 treated 
above. The plots of the effective central charge given in Fig. |lj| again give us strong 
reasons to believe that the a —>■ oo limit represents the non-trivial critical fixed point 
predicted by perturbation theory. Also the effective thermal and magnetic exponents of 



Figs. |TT|.b and |12|.b respectively exhibit a structure similar to the q = 3 case, but the 



critical exponents at a — > oo are of course distinct. 
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Figure 11: Thermal exponent for x t for three coupled g-state Potts models with q = 3 
and q = 4 respectively. 
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Figure 12: Magnetic exponents x$ and x$ for three coupled g-state Potts models 
with q = 3 and q = 4 respectively. 
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Figure 13: Central charge for three coupled 4-state Potts models. 
5.3.4 The a — ► oo limit for general g G [2,4] 

After having numerically identified the non-trivial critical fixed point for three coupled 
models, that was predicted by the perturbative RG, we would like to accurately determine 
the associated critical exponents as a function of q. This is done by diagonalising the 
transfer matrix directly at a = oo, using alg4. 



In Table |ll] we show finite-size estimates of the central charge, obtained from three- 
point fits of the form ( [4.81 ). For all q G [2,4] these estimates seems to converge as 
L — > oo. Of course, an explicit extrapolation to the infinite system limit from just three 
data points would easily turn out to be rather subjective. It should however be noted 
that the convergence is in most cases monotonic, and the last estimate therefore provides 
a quite accurate upper or lower bound for the extrapolated quantity. 

The results for c are in very good agreement with the perturbative RG in the range 
q G [2, 3]. For larger values of q the perturbation theory is expected to break down, and 
so the numerical values are likely to be more reliable. 



The corresponding values of the thermal scaling dimension are given in Table [12[ In 
this case the comparison with the perturbative results is more tricky. We recall that x t 
is obtained by fitting the first gap in the even sector of the transfer matrix to Eq. (|4.9|) . 
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Q 


c(4,8) 


c(6, 10) 


c(8,12) 


2.00 


1.489114 


1.497443 


1.500008 


2.25 


1.747538 


1.757162 


1.760185 


2.50 


1.976384 


1.986719 


1.989905 


2.75 


2.179933 


2.190166 


2.193005 


3.00 


2.361502 


2.370594 


2.372345 


3.25 


2.523740 


2.530453 


2.530147 


3.50 


2.668818 


2.671752 


2.668217 


3.75 


2.798557 


2.796195 


2.788092 


4.00 


2.914511 


2.905265 


2.891139 



Table 11: Results for the central charge of three coupled models at a = oo and for varying 
values of q. 



1 


x t (6,8) 


x t (8,10) 


^(10,12) 


2.00 


1.257900 


1.248706 


1.242186 


2.25 


1.263919 


1.255942 


1.250623 


2.50 


1.268895 


1.262901 


1.259341 


2.75 


1.272294 


1.268558 


1.266905 


3.00 


1.273618 


1.271938 


1.271901 


3.25 


1.272423 


1.272143 


1.272966 


3.50 


1.286448 


1.268344 


1.268400 


3.75 


1.261129 


1.260138 


1.258753 


4.00 


1.250663 


1.247052 


1.242082 



Table 12: Results for the thermal scaling dimension of three coupled models at a = oo 
using alg4. 
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Q 


c(3,5) 


c(4,6) 


x t (4,5) 


x t (5,6) 


2.00 


1.5040 


1.4906 


1.2692 


1.2600 


2.25 


1.7653 


1.7490 


1.2798 


1.2691 


2.50 


1.9965 


1.9775 


1.2890 


1.2774 


2.75 


2.2022 


2.1801 


1.2964 


1.2842 


3.00 


2.3856 


2.3601 


1.3018 


1.2890 


3.25 


2.5495 


2.5198 


1.3051 


1.2913 


3.50 


2.6962 


2.6615 


1.3062 


1.2907 


3.75 


2.8274 


2.7869 


1.3049 


1.2868 


4.00 


2.9449 


2.8975 


1.3011 


1.2796 



Table 13: Results for three models at a — > oo using alg3. 



This is expected to correspond to the scaling dimension of the most relevant symmetric 
energy operator, E\ + 62 + £3. However, A ei+£2+£3 has only been determined to order 
0(e 3 ) [see Eq. 03.11 )1, and so the discrepancy between the values of Table and of 
Table |l] for large q should hardly come as a surprise. On the other hand, for q slightly 
larger than two numerics is supposed to be hampered by large logarithmic corrections 
due to the near-marginality of the energy operator ||40|| . In this range we thus expect the 
perturbative results to be the more reliable. 



In Table |T~3| we show some more results for c and x t , but this time obtained from alg3. 
The usefulness of these data is twofold. First, they clearly demonstrate the superiority 
of alg4 as regards the accuracy and the convergence properties (monotonicity) of the 



finite-size estimates. Second, the agreement with Table 11-12 serves as a check of the 



modular invariance of the a = 00 critical point, as explained in Sec. 4.7. 



We now turn our attention to the magnetic exponents x^ mag ' 1 obtained by imposing 
twisted boundary conditions on iV ma g = 1,2 or 3 of the N = 3 coupled Potts models. 
The numerical results are presented in Table |T4]-|16[ Comparing with the analytic results 
of Table [T]-^ we find, as expected, that these scaling exponents must be those of the 
first three moments of the magnetisation. Excluding the case of q = 2 where the remarks 
made above hold true, we find a very good agreement between Xjj and A CT1 , and between 
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q 


xS } (4,6) 


4 } (6,8) 


4? (8, 10) 


2.00 


0.124815 


0.125054 


0.125112 


2.25 


0.127712 


0.127847 


0.127851 


2.50 


0.129955 


0.129915 


0.129810 


2.75 


0.131635 


0.131334 


0.131050 


3.00 


0.132825 


0.132170 


0.131623 


3.25 


0.133588 


0.132480 


0.131577 


3.50 


0.133977 


0.132318 


0.130962 


3.75 


0.134040 


0.131734 


0.129831 


4.00 


0.133816 


0.130778 


0.128237 



Table 14: x^ for three coupled Potts models at a = oo. 



(2) 

x H and Ao- l0 - 2 in the range q < 3. For larger q the perturbative RG breaks down with a 
vengeance, as witnessed by its predictions of negative scaling dimension. 

Finally, the agreement between x H and A aif72Cr;i , especially in the range 2.5 < q < 3.5, 
is quite striking, and is the first validation of the second order perturbative computation 
proposed previously by one of us |2"0"|. 



5.3.5 Higher exponents in the even sector for N = 3, q = 3 

In an attempt to numerically identify more of the operators predicted by perturbation 
theory, we have looked at the scaling dimensions extracted from the finite-size scaling of 
higher gaps in the even sector of the transfer matrix. These computations are extremely 
time-consuming, since, in order to study the first k gaps, we need to iterate and orthog- 
onalise k + 1 vectors, cfr. Eq. ( }4.3|) . We have therefore focused exclusively on the case of 
three coupled 3-state Potts models. 

The primary operators we expect to find in the even sector must all be energetic, 
since the Z q symmetry associated with the permutation of the Potts spin labels has been 
factored out. Furthermore, since alg4 by construction treats all three layers symmetri- 
cally, we only expect to find such operators that are symmetric under any permutation 
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q 


4* (4,6) 

STL V " / 


4? (6, 8) 


4? (8, 10) 


2.00 


0.277424 


0.276866 




2.25 


0.287411 


0.287143 


0.286780 


2.50 


0.296026 


0.295985 


0.295839 


2.75 


0.303408 


0.303470 


0.303445 


3.00 


0.309671 


0.309661 


0.309600 


3.25 


0.314915 


0.314615 


0.314310 


3.50 


0.319226 


0.318393 


0.317599 


3.75 


0.322685 


0.321061 


0.319512 


4.00 


0.325367 


0.322693 


0.320119 



Table 15: x H for three coupled Potts models at a = oo. 



Q 


of (4, 6) 


4? (6, 8) 


4^(8,10) 


2.00 


0.475424 


0.471563 




2.25 


0.503049 


0.500930 




2.50 


0.529688 


0.529874 




2.75 


0.555511 


0.558541 




3.00 


0.580638 


0.587025 


0.592223 


3.25 


0.605162 


0.615388 




3.50 


0.629150 


0.643668 




3.75 


0.652656 


0.671885 




4.00 


0.675719 


0.700047 


0.721371 



Table 16: x H for three coupled Potts models at a = oo. 
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Gap 


x(4) 


x(6) 


x(8) 


x(10) 


x(12) 


Extrapolation 


Operator 


1 


1.694 


1.471 


1.385 


1.344 


1.322 


1.27 


£i + e 2 + £3 


2 


2.603 


2.380 


2.272 


2.219 


2.148 


« 2.0 


T = L_i7 


3 


2.922 


2.775 


2.398 


2.225 


2.148 


« 2.0 


T = L_ x l 


4 


— 


3.104 


2.398 


2.225 


2.174 


« 2.1 


E\E2 + ^2^3 + ^3^1 


5 




3.104 


2.685 


2.633 


2.599 


« 2.3 


L_i(ei +£2 + ^3) 


6 




3.626 


3.279 


2.785 


2.598 


« 2.3 


L-iOn +£2 + ^3) 


7 




3.991 


3.279 


2.785 


2.598 


« 2.4 


£l£2£3 


8 




3.991 


3.289 


3.146 


3.066 


« 3.0 





Table 17: Scaling dimensions extracted from the higher gaps in the even sector as well 
as their identification with physical operators. 



of the layer indices. According to Sec. ^| the permissible operators are thus 

J, E\ + 62 + £3; ^1^2 + ^2^3 + ^3^1; £l£2£3- 



(5.9) 



If the (presently unknown) conformal field theory of the system can be assumed to have 
the usual structure |42| we can, in addition to these primaries, in general expect to 
find descendent operators with their appropriate degeneracies, reflecting the null vector 
structure of the Verma module. The identity operator I can be thought of as the primary 
corresponding to the zeroth gap, thus having scaling dimension Aj = by definition. 

In Table O we show the results for the first eight gaps, obtained by employing alg4. 
Since we do not wish to make any a priori assumptions on the existence of the stress 
tensor we report only the unadorned one-point fits ( f4.5|) , for strip widths up to L = 12. 
For L = 4 we were only able to iterate four orthogonal vectors, and consequently there 
are some empty entries in the table. 

The extrapolation to the L — > 00 limit becomes increasingly important for the higher 
gaps, but unfortunately also increasingly difficult, since the convergence is slower and 
there are fewer data points. The reported values of the extrapolated scaling dimensions 
are therefore only indicative. 
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Let us recall from Sec. |3| that perturbation theory predicts 

A £1+£2+£3 ~ 1.44, A £l£2+£2£3+£3£l ~ 2.08, A £l£2£3 ~ 2.28. (5.10) 

These values cannot be believed to be very precise, since at q = 3 we expect the pertur- 
bative expansion to be at the verge of breaking down. Nevertheless, they are supposed 
to be of the right order of magnitude, and we thus venture to make the identification 
shown in Table [17]. 

Note that for L > 8 the 3rd and the 4th gap, and also the 6th and the 7th, are exactly 
degenerate with full 16-digit machine precision. One would therefore expect them to be 
connected by a simple symmetry, as reflected by our conjectured identification. However, 
for L = 12 a crossover takes place so that henceforth it is the second and the third gap that 
are degenerate rather than the third and the fourth. This could have been anticipated 
by noticing that for L < 12 gaps 2-4 have quite different finite-size dependence. A 
similar crossover is expected to take place amongst gaps 5-7 for L = 14. Accordingly our 
operator identifications pertain to the expected arrangement of the gaps in the L — > oo 
limit. 

Apart from the primaries ( |5.9|) we are able to identify the holomorphic and antiholo- 
morphic components of the stress tensor, both with scaling dimension two. This provides 
further evidence for the criticality of the system, and is completely independent of the 
methods used this far. In addition we observe the level-one descendents of the operator 
£i + £2 + £3 with the expected scaling dimension. 

Unfortunately we have not been able to iterate enough vectors to see descendents at 
level two. This is a pity, since the degeneracy of these scaling dimensions would enable 
us to make predictions about the null vector structure of the Verma module. 

5.3.6 Monte Carlo Results 

Showing the existence of scaling operators, whose correlators obey scaling laws, would 
confirm that the a — > 00 limit point on the self-dual line for three coupled models can 
be identified with a second order phase transition. To do so, we performed Monte Carlo 
(MC) simulations with the simplified Boltzmann weights (|2.36|) pertaining to this point. 
A direct measure of correlators for the physically significant operators (i.e. those which 
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are diagonal in the perturbation scheme) clearly shows that the model is indeed critical, 
exhibiting scaling laws with associated exponents sufficiently close to those obtained 
either by perturbative CFT or transfer matrix techniques. Since the critical exponents 
can be measured with much more accuracy using transfer matrices techniques, MC was 
used solely to establish the criticality of the fixed point model. 

We chose to measure explicitly the correlation functions for the spin and energy oper- 
ators. Whilst the definition of energy correlators is straightforward, the Z q symmetry of 
spin variables has to be taken into account, and one should rather consider the operators 

Ej = <Xj — W{ mod q, (5-H) 

where Wl is the instantaneous average value of the spin o on lattice i. The physically 
meaningful operators, for which we computed correlators, are thus the following: 

e s = Sx + 62 + e 3 (5.12) 

e A = e 1 -e 2 (5.13) 

£ 5 = Sx + S 2 + S 3 (5.14) 

£4 = Si-E 2 . (5.15) 

These operators have well-behaved distributions and are the ones obeying scaling laws. 
We performed simulations on a 160 x 160 lattice with periodic boundary conditions 
in both directions. Sixteen runs of 10,000 sweeps each were made, with a preliminary 
thermalisation stage of 500 times the autocorrelation time (which is around 25 updates 
for our lattice size). We also checked that the values of critical exponents were rather 
close to those predicted using this lattice size. A larger lattice would evidently provide 
a better description of critical behaviour, but computation time grows rapidly with the 
number of spins and, as we said, a precise measurement of critical exponents is not what 
was sought here. 

Simulations were performed using a modified version of the Wolff one-cluster algo- 
rithm [H3J. Let us consider the model given by Eq. (2.16), which we recall: 



Hij = a (5 a . >crj + S n >T . + S Vi >v . ) + b(5 at >a , 5 n >T . + 5 Ut >a . S Vi m + 6 n >Tj 5 m >Vj ) 

+ c <W.<WAi,%- ( 5 - 16 ) 
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We first choose randomly one of the three models, say the one with the indices a (the 
model is invariant under permutation of any pair of layers). Then, we consider the model 
defined by 

Hij = (a + b (5 Ti>Tj +6 Vurh ) + c 5^5^)5^,^ + const. (5.17) 
= a'8 a . )(r . + const. (5.18) 

Next we update the lattice with the usual Wolff algorithm but with the local a' = 

a + b (S Ti}Tj +S VitVj ) +c &Ti,Tj&mMi • Repeated three times, this operation defines one update 
of the system. 

The measurement procedure is quite straightforward: For a given operator O, all cor- 
relation functions G(\x—y\) = (0(x)0(y)) , which are manifestly translationally invariant 
due to the choice of boundary conditions, are computed for < \x — y\ < 80. Compiling 
results, one finds that the region where scaling laws exist is limited to \x — y\ < 20 for 
all operators, except for the symmetric energy, whose scaling law behaviour is observed 
only for \x — y\ < 8, this being due to the fact that the critical exponent for this operator 
is much larger than for the others. This is quite deceiving, but not unexpected. It is 
clear that for critical correlations to have some sense, they must be between points x and 
y such that \x — y\ <C L. This, combined with the fact that our statistics were rather 
low, explains why we had to restrict our attention to such narrow regions. Doing so, one 
observes nice scaling properties with associated exponents close to the ones computed 



by other analytical and numerical methods. Figures [LJ and [15] present the Log-Log plot 
of correlation functions along with statistical error bars. Linear fits of these graphs give 
2Ao, and thus leads to the following values for the critical exponents: 

A es = 1.27 ±0.13 (5.19) 

A £A = 0.63 ±0.04 (5.20) 

A Ss = 0.13 ±0.03 (5.21) 

A Ea = 0.10 ±0.03. (5.22) 

Altough not precise, the MC data clearly establish the existence of scaling operators. 
This confirms that the model under study is critical. 
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Figure 14: (a) Monte Carlo results for \og(s a(x)s a(u)) as a function of log \x — y\. (b) 
log (es(x) eg (y)) as a function of log \x — y\. 
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Figure 15: (a) Monte Carlo results for log(XU(:r)£yi(2/)) as a function of log |x — y\. (b) 
log(Es(x)Eg(y)) as a function of log \x — y\. 
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6 Concluding remarks 



Besides the long-reaching goal of describing exactly disordered systems, the extensive 
study presented here is, we hope, interesting in many aspects. First, it introduces a 
new variant of the Potts model transfer matrix method which proves to be the most 
accurate and efficient up to now. Second, we have identified non-trivial fixed points both 
numerically and analytically and have presented evidence that these two types of fixed 
points are indeed the same and correspond to non-trivial critical models. The universality 
classes of these critical models are new, and several critical exponents were computed for 
the first time. 

The next step in our long-term project is to solve the critical model, which was proven 
to be the end-point of the self-dual line. At this point, as we have shown, the Hamiltonian 
simplifies greatly, and it might be possible to map the system to a loop model, whose 



continuum limit is a Liouville field theory [36]. We believe that the results presented 
here will be of the outmost importance when searching for such a loop formulation. 
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